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ABSTRACT 

Forward and aft acoustic propagation and radiation from a ducted fan is modelled 
using a finite element discretization of the acoustic field equations. The fan noise source is 
introduced as equivalent body forces representing distributed blade loading. The flow in 
and around the nacelle is assumed to be nonuniform, reflecting the effects of forward flight 
and flow into the inlet. Refraction due to the fan exit jet shear layer is not represented. 
Acoustic treatment on the inlet and exhaust duct surfaces provides a mechanism for 
attenuation. In a region enclosing the fan a pressure formulation is used with the assumption 
of locally uniform flow. Away from the fan a velocity potential formulation is used and the 
flow is assumed nonuniform but irrotational. A procedure is developed for matching the two 
regions by making use of local duct modal amplitudes as transition state variables and 
determining the amplitudes by enforcing natural boundary conditions at the interface 
between adjacent regions in which pressure and velocity potential are used. Simple models 
of rotor alone and rotor/exit guide vane generated noise are used to demonstrate the 
calculation of the radiated acoustic field and to show the effect of acoustic treatment. The 
model has been used to asses the success of four techniques for acoustic lining optimization 
in reducing far field noise. 

INTRODUCTION 

Figure 1 shows in idealized form a rotor/exit guide vane configuration imbedded in 
a nacelle with a centerbody or a core engine. The rotor represents the fan in a high bypass 
turbofan engine or a ducted propeller. The exit guide vanes provide a source of interaction 
noise. In the numerical examples considered in this investigation the number of blades and 
exit guide vanes is characteristic of a ducted propeller. The rotor/exit guide vane source 
generates noise which is propagated through the inlet and exhaust ducts and is radiated to 
the far field. Nonuniform steady flow exists in and around the nacelle due to inflow, outflow, 
and forward flight effects. It is required to predict the far field radiated noise at harmonics 
of the blade passage frequency. 

In a previous investigation [1] ducted fan noise was studied with the assumption that 
the mean flow in and around the inlet could be assumed to be uniform. This allowed the 
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acoustic field equations to be simplified to a convected wave equation in the acoustic 
pressure. The fan noise source was introduced as equivalent body forces representing 
distributed blade loading. In the investigation reported here the flow in and around the 
nacelle is assumed to be nonuniform, reflecting the effects of forward flight and flow into 
the inlet. Refraction due to the fan exit jet shear layer is not accounted for in the current 
model. Nonuniform mean flow eliminates the possibility of using the convected wave 
equation as in Reference [1]. Nonuniform flow effects have been included in studies of 
forward radiated noise from turbofan inlets [2,3]. In this case the mean flow is assumed to 
be irrotational, as is the acoustic perturbation, and the introduction of a velocity potential 
simplifies the field equations. The fan noise source is introduced as a boundary condition 
on the fan face defining the amplitude of incident and reflected duct modes. 

In order to ensure an efficient numerical model it is desirable to represent the 
acoustic field in terms of a single state variable. In the previous studies this was done in 
terms of the acoustic pressure or the velocity potential. In the extension discussed here the 
pressure formulation is not suitable because the mean flow is nonuniform in most of the 
region of propagation and radiation. The velocity potential formulation is not suitable for 
introducing the fan noise source as an equivalent body force distribution. For these reasons 
a mixed method is introduced. In the fan region a pressure formulation is used and away 
from the fan a velocity potential formulation is used. When pressure is used as the state 
variable in a region near the fan where it is assumed possible to consider a locally uniform 
flow, the fan noise source, either rotor alone noise, or interaction noise, can be introduced 
as equivalent body forces. In regions away from the fan where the mean flow is nonuniform 
but irrotational, a velocity potential description has proven to be appropriate. A procedure 
is developed for matching the two regions by making use of local duct modal amplitudes as 
transition state variables and determining the amplitudes by enforcing natural boundary 
conditions at the interface between adjacent regions in which pressure and velocity potential 
are used. 

An additional feature introduced in the present study is the provision for modelling 
an acoustic lining on the surfaces of the fan inlet and exhaust ducts. The lining is assumed 
to be point reacting and to have a frequency dependent impedance or admittance. The 
lining representation is consistent with the physical characteristics of current acoustical 
materials. 

The details of the introduction of the noise source have been discussed in [1]. In this 
paper the extension of the formulation of the finite element model required to admit 
nonuniform mean flow and acoustic treatment is explained. Results are presented for the 
far field radiated acoustic field for examples of rotor alone noise and for rotor/EGV 
interaction noise to demonstrate the cut off of subsonic tip speed rotor alone noise, the 
propagation of interaction tones, and the influence of acoustic treatment on the radiated 
field. 

PROBLEM FORMULATION 

In previous studies acoustic radiation from fan sources embedded in a shroud or duct 
has been modeled using two different formulations. In the case of forward radiated noise 
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from turbofan nacelles [2,3] the model was based on the assumption that the mean flow in 
and around the nacelle is irrotational and that the acoustic perturbation is also irrotational. 
This makes it possible to introduce mean flow and acoustic perturbation velocity potentials. 
The governing field equations are the linearized acoustic continuity equation and a 
linearized acoustic Bernoulli equation, written in terms of the acoustic potential and the 
acoustic pressure (or density). As shown in [2,3], the acoustic potential is the solution of 


^tV-(p r V4> + pV4,,) = 0 
ot 


( 1 ) 


and 


P 


- -%-^ + V<tyV4>) 


X dt 


( 2 ) 


where <f> is the acoustic potential, <£ r is the local mean flow (reference) potential, p is the 
acoustic density, p r is the local mean flow density, and c r is the local speed of sound in the 
mean flow. All quantities are nondimensional with respect to the density in the far field, 
Pa, , the speed of sound in the far field, Co, , and a reference length, R, which is the fan 
radius. The acoustic potential is nondimensional with respect to Ca,R , and the acoustic 

pressure with respect to p.e. . Time is scaled with R/c* . The fan or EGV source is input 
by specifying complex amplitudes of duct modes at a boundary of the computational domain 
designated as the source plane. In the acoustic potential formulation the option to describe 
the source in terms of equivalent volumetric forces within the computational domain does 
not appear to be available because of the use of the acoustic Bernoulli equation. The 
acoustic field equations (1) and (2) form the basis of the finite element models for acoustic 
radiation from turbofan inlets including the effect of forward flight developed in [2,3]. 
Comparisons of results with flight test data are described in [4]. 

Acoustic radiation from unshrouded propellers in a free field and in a wind tunnel 
environment has been investigated with the assumption that the mean flow field, 
representing the forward flight effect or the flow in the wind tunnel, is nearly uniform [5-8]. 
In this case the governing acoustic field equation is the convected wave equation in terms 
of the acoustic pressure 


V-(Vp- 


— dp i nf\ 'l M r 1 &P 

c *dx p/) c *a*t c * dt 2 


= 0 


(3) 
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The fan or EGV source is introduced by an equivalent distribution of body forces per unit 
mass / (referred to herein as volumetric force). Equation (3) is in nondimensional form 
with the body force nondimensionalized by cJR . The volumetric forces have been derived 
from a simplified lifting line theory in Reference [1], The appearance of c T is due to the 
use of c M in the definition of dimensionless variables. MJc r is the local Mach number M f . 
The rotor alone model considers the steady blade loading rotating with the blades, while the 
EGV noise considers stationary blades with unsteady loading produced by the sweeping of 
rotor wakes past the blades. While simplified in the details, this approach contains the 
essential physics of two principal types of turbomachinery noise. 

The extension of the propeller acoustic radiation model to shrouded and ducted fans 
was done with the assumption that in this case also the entire flow field, both internal and 
external to the nacelle, is uniform [1]. This is probably acceptable for a shrouded propeller 
and for a very high bypass ratio ducted fan. For a more typical turbofan nacelle the 
assumption is clearly not satisfied. The physical appeal of the source model motivates the 
extension of the concept to include nonuniform mean flow. The virtue of the representation 
of the source is that it provides a direct link between blade loading and acoustic source 
strength and distribution. It is not required to specify the acoustic duct modal amplitudes 
produced by the source. 

The approach used is a mixed method in which the region containing the noise 
source is a pressure formulation based on the convected wave equation (3), and away from 
the source the acoustic potential formulation of equations (1) and (2) provides the field 
equations. The only assumption required is that within the narrow source region it is 
sufficiently accurate to treat the mean flow as uniform. Figure.,2 is a sketch of a fan or EGV 
blade row embedded in a nacelle, with flow moving from right to left. A finite element mesh 
based on quadrilateral elements in the axially symmetric geometry is shown superimposed 
on the nacelle fan inlet and fan exhaust ducts. Of particular interest in this discussion is a 
region near the blade row consisting of four columns of elements in which the field equation 
is based on pressure. The shaded column contains the acoustic volumetric forces describing 
the blade loading. In this narrow region the flow is taken as axially directed and uniform. 
Outside this region the potential formulation is used. The details of the finite element 
formulation in the pressure region are given in [5,6] and the corresponding development for 
the potential formulation is found in [2,3]. Reference [1] should be consulted for the source 
model. 

The main issue to be resolved is the matching between the two regions. In order to 
see how this can be done it is necessary to reiterate the finite element procedure for the 
weak formulation of the two types of problems. In both cases a Galerkin weighted residual 
method is used. In the case of the convected wave equation a solution for the pressure field 

p(xs)e' (r, ' x ~ me) is sought such that the weighted residual statement 
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(4) 


fffw{V<Vp- 


Mi 


.2 dx 




frl-n A-' 


+ = 0 
c r 


for every member W(x,r)e‘ (r ' ,t * md) of a complete set of functions. In equation (4) the 
harmonic time dependence of the source is explicitly represented by the nondimensional 
frequency rj r defined as uR/c n . In the form shown the operations implied require that 
W(x,r) be piecewise continuous and that p(x,r) have a continuous first derivative, 
requirements which would lead to seeking solutions in a very restrictive class of functions. 
A weaker solution, one which admits the possibility of a solution in a less restrictive class 
of functions is obtained by integration by parts to yield the weighted residuals statement that 

solutions p(xf)e from the class of continuous functions are sought which satisfy 




- ff W(Vp- ——0'ndS = 0 

i J r 1 dx 


(5) 


for every member W(xf)e ,(n,t * md) of a set of continuous functions. The volume integral 
is over the domain in which the pressure formulation is used, and the surface integral is over 

the boundary of this domain. The unit normal n is directed out of the domain. The 
boundary integral introduced plays a significant role in the matching procedure. 

In the case of the potential region the weak formulation seeks solutions 

4 >(x,r)e ,(v ~* ,e) in the class of continuous functions which satisfy the weighted residual 
equation 


/// 1 W<p,V4. V<t>,p) -iri^p ) dV 

V 

-ff W(p r V<f> + V4> r p ) mdS = 0 


in which p(x,r)e lt ' n '*~ m9) is defined by equation (2), for every member W{x,r)e l(A,t * m&) 
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of a set of continuous functions. The boundary integral is important in the matching 
procedure. 

The boundary integrals represent natural boundary conditions which must be imposed 
on the boundary of the domain. The use of the boundary integrals on the boundaries of the 
problem is discussed in [2]. These boundaries of the computational domain are shown in 
Figure 3. The far field boundary is at a large distance from the nacelle and is a non 
reflecting surface on which a radiation condition is applied [2,3]. This surface is the outer 
boundary of wave envelope elements which allow a transition from a fine mesh near the 
nacelle to a very coarse mesh in the far field. Most of the nacelle and centerbody surfaces 
are rigid, where the normal component of acoustic particle velocity vanishes. In the pressure 

region this corresponds to — = 0 . With the additional condition that the flow is taken to 

dn 

be axially directed in the pressure region, it follows that i-n = 0 , so that the boundary 

integral vanishes. In the velocity potential region — = 0 . In addition, the flow tangency 

dn 

condition requires that V<j> r -n = 0 . The boundary integral vanishes in this case as well. A 
portion of the nacelle and centerbody is acoustically treated. On these surfaces an 
impedance relation is specified, as discussed later. 

In the mixed formulation there are two domains, one in which pressure is the 
dependent variable, and one in which velocity potential is the dependent variable. The 
boundary integrals at the interface of the two domains represent the natural boundary 
conditions which one domain imposes on the other. For the pressure region this integral is 



(7) 


The positive sign applies if the pressure region boundary has its normal in the positive x 
direction. The velocity potential region boundary integral is 


*v - *)dS 


dx 


( 8 ) 


The sign choice is the same as noted for the pressure case. These integrals will form 
contributions to the "stiffness" matrices obtained in the finite element formulation for the 
elements on the boundary between the regions. 
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MATCHING OF SOLUTIONS 


The matching at the interface between the two domains is accomplished by using the 
connection between pressure and velocity potential provided by combining equation (2) and 
the linearized equation of state p = c r 2 p to yield 


P = -P r 0'h r 4> + M r -|£) = -? T c r (ir\fi + M f Q) (9) 

The nondimensional frequency r\ f and Mach number M { are defined relative to the local 
conditions at the interface, taken as those at the fan source. 

At the interface the velocity potential is written as an expansion in terms of the local 
acoustic modes for the duct, 


M 


' kli * + b$~ i (r)e '* x,X ] 


i-1 


( 10 ) 


where the duct eigenfunctions \|rf(r) are computed using a finite element formulation at 
the interface cross section. 

The duct eigenfunctions .(r) are solutions of the Bessel equation and boundary 
conditions on the duct wall 


1^0 

r dr dr r 2 

= 0 r = „ (11) 

dr 

=0 r - 1 

dr 


a is the nondimensional inner radius of the duct. The duct eigenfunctions are the same for 
upstream propagation and downstream propagation so that ilqCr) = i|f,"(r) = ijr^r) . The 
axial wave numbers k* are calculated from the duct eigenvalues k, according to 
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( 12 ) 


1/ (1 -«/) 


-M,± 


N 


1-(1 - M ^) (— ) 2 
1/ 


where the mach number M } and the nondimensional frequency ri r are local values at the 
interface cross section. The eigenvalue problem posed by equations (11) and (12) is solved 
utilizing a finite element formulation on a one dimensional grid of quadratic elements which 
matches the grid used in the duct interior. This very robust routine produces eigenvalues and 
eigenfunctions of high accuracy which are completely consistent with the interior grid. 
Reference [9] provides some details of the procedure. Equation (9) then provides an 
expansion for the pressure at the interface in terms of the velocity potential modal 
amplitudes 


M 


kt 


P = -JPrWE +b^(r)(l-M f -^)e ,<v ] 


K 


-ifc'x 


1=1 


(13) 


Equations (10) and (13) also provide expansions for terms in the integrands of 
equations (7) and (8) 




M 


i = 1 


S ti. 


lc jfc 
*1/ 


-**> 


(14) 


M? 54, . M 


M 

-*Pr 1 

i-1 


(1 -M*X — •) 


i|f t .(r)e 




+ b i 


. k~ 

(l-M/)(-l)*M. 
1/ J 




( 15 ) 


The expansions are written in vector-matrix form as 
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(16) 


<t> = I '!'<('■) Je] 


n 

\b\ 



(17) 


2 ( > 

(1 ~“^)& = "*lwiWK Y i e O 

c r l J 


(18) 


(l-*)* 
c y ex 



<t> = -i[$i(r) | !•(/■)][ a le 



(19) 


The row matrix [t|r(r)] = [\jr f (r) | ijr.(r)] has 2M columns constructed from the continuous 
eigenfunctions generated by the eigenvalue problem of equation (11). It is partitioned with 
two blocks of M columns of the eigenfunctions retained in' the expansion. The diagonal 

square matrix [e] of size 2M x 2M has elements e» ,* = 1 M and e~ ,i = M+12.M , 
where 






( 20 ) 


The diagonal square matrix [a] has elements a. I ,i - 1 M , and ,i = M+l£M , 
where 


K 


a„- = p,iy(l 


( 21 ) 
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The diagonal square matrix [P] has elements P^ ,i=lM , and ,i=M + l£M , where 



-ip 



( 22 ) 


Finally, the diagonal square matrix [y] has elements y~ ,i=lM , and y~ ,i=M+l£M , 
where 


it* fc 1 

vj = -ip,C,T|/(l -J#yr-5)(— ) 

1 / 1 / 


(23) 


The vector 



is partitioned with the complex amplitudes a. 


i=l,M for right running 


acoustic modes and b i , i=l,M for left running modes. 

The finite element implementation of the Galerkin Method of Weighted Residuals 
(MWR) is characterized by interpolation within subdomains of the computational domain 
(elements) based on values of the field variable, pressure or velocity potential, at nodes of 
the elements. In the work reported here the elements are isoparametric quadrilaterals with 
eight nodes. This type of discretization provides continuity of the field variable at element 
boundaries, but does not produce solutions with continuous derivatives. This is consistent 
with the weak formulation. Each element contributes an element "stiffness matrix" to the set 
of linear equations whose solution yields the nodal values of the field variable. 

Figure 4 shows the interface between a region of pressure and a region of velocity 
potential. The elements in the two regions which are on the boundary are the key to the 
matching of the solutions in the two regions. In Figure 4 the pressure elements have their 
right boundary on the interface and the potential elements have their left boundary on the 
interface. The continuity of the solution across the interface is accomplished by using the 
fact that both pressure and velocity potential at the interface can be expanded in terms of 
the duct acoustic modes appropriate to the geometry and flow conditions at the interface 
as described in equations (16) and (17). These finite eigenfunction expansions contain M 
duct modes in each direction. The value of M is chosen to include all propagating modes 
plus perhaps three to five cutoff modes to assure radial resolution of the pressure field. 

The velocity potential element stiffness matrices on the boundary are transformed 
by using the discrete eigenvectors to form a transformation matrix such that 


10 



a 


(24) 


W - M 


where {<j>} , the vector of nodal values of the velocity potential for an element, are replaced 
by the modal amplitudes a ( , b. and the nodal values of <J> at the interior nodes <j>. in the 
element (nodes not on the boundary). A similar transformation is made for pressure 
elements 




« = 


(25) 


The elements of [Y x ] and [Y 2 ] are constructed by using the discrete eigenvectors, sampled 
at the boundary nodes for the element, and by using the eigenvector expansions of equations 
(16) and (17). They are "modal matrices" which serve as transformations from nodal values 
of the field variable in the interior of the element and modal amplitudes to values of the 
field variable at all of the nodes. A similar interpretation can be given to an eigenfunction 
expansion of the weighting functions. The weighting function evaluated at the element nodes { W) 
is obtained in velocity potential elements as 


m - ra 


a 

b 


(26) 


and in pressure elements as 


m = m 


a 

b 


( 27 ) 


The transformed element stiffness matrices are for pressure 
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M - 


(28) 


and for the velocity potential 

M ■ ['•M'fcl*.] 


(29) 


The right boundary of the pressure elements and the left boundary of the potential elements 
at the interface are now discretized in terms of the modal amplitudes ,b l . The assembly 
process is carried out on the basis that a. and b t on the pressure and potential elements are 
the same. 

The element stiffness matrices on the boundary must be augmented by the addition 
of the boundary integrals. These integrals can be conveniently evaluated directly in terms 
of the modal amplitudes and interior nodal values. This is accomplished by using the 
eigenfunction expansions for the integrands given by equations (18) and (19). For example, 
consider the integral of equation (7) evaluated on the right boundary. It can be written 



The integral of equation (8) evaluated on the left boundary is similarly written 

/„ - +(-o//w[*(r)I«]^J<iS (31) 


The eigenfunctions ^(r) , treated in the expansions as continuous functions of r, are known 
only in terms of their discrete values (eigenvectors) at the nodal points on the interface in 
the finite element eigenvalue problem described by equation (11). The row matrix of 
continuous eigenfunctions is obtained by interpolation of the corresponding discrete modal 
matrix [Y] = [Y. | Y ( .j of discrete eigenvectors according to 

[*<f>] - [NJT] 
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( 32 ) 



where [N b ] is the element interpolation matrix on the right boundary of a pressure element 
or left boundary of a potential element. In standard finite element procedure, [ N b ] 
produces a continuous version of a function on the boundary in terms of discrete values of 
the function at nodes on the boundary. The continuous weighting functions can also be 
obtained from discrete values at the nodes { W] . This is accomplished using the same 
interpolation matrix 

W - [tf i ](W'} (33) 


For pressure element weighting functions 


{W) = [TJP 


fl 


and for potential weighting functions 


(34) 


iW) -[T 


fl 


(35) 


The boundary integrals / and / v can be written as 




(36) 


I, - -'•[5P] r //[W 4 ] T [W s )iS[TIa|} 


(37) 


The element integral contributions are added to the element stiffness matrices prior to 
assembly. The addition of the boundary 
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integral contributions to the stiffness matrices for elements on the boundary can be 
simplified if additional "modal matrices" [*P 3 ] and [Y 4 ] , constructed by using the discrete 
eigenvectors, sampled at the boundary nodes for the element, and by using the eigenvector 
expansions of equations (18) and (19), are introduced to relate internal nodal values of the 
field variables and modal amplitudes on the boundary to values of the field variables at the 
element nodes. Equations (36) and (37) can be rewritten as 




V 

a 


(38) 




a 

b\ 

H>, 


(39) 


The interpolation matrix [ N ] is the element interpolation matrix which produces a 
continuous function within an element in terms of the nodal values of the function. The 
advantage of this formulation is that the boundary contribution is the same dimension as the 
transformed stiffness matrix to which it is appended. 

ACOUSTIC TREATMENT 

No previously published studies of fan noise radiation have addressed the placement 
of acoustic treatment on the surfaces of the fan inlet and fan exhaust ducts. In the 
formulation described here provision has been made for acoustic treatment in the region 
in which the acoustic field is described in terms of the velocity potential. This excludes only 
a very small region in which the fan noise generation process is described in terms of a 
pressure formulation. 

The boundary integral of equation (6) is the mechanism by which the boundary 
condition imposed by locally reacting acoustic treatment is introduced. On surfaces on which 
acoustic treatment is present the normal component of mean flow velocity vanishes and the 
lining boundary integral simplifies to 



( 40 ) 
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where v n is the acoustic particle velocity directed normally into the acoustic treatment. The 
acoustic treatment is described in terms of the impedance relationship 


P 

vn 



(41) 


p is the acoustic pressure and vn is the normal component of lining velocity at the wall. 
The impedance z is a prescribed function of frequency and is nondimensional with respect 
to • A is defined as the nondimensional acoustic admittance. The relation between the 
fluid particle velocity at the wall and the wall velocity is one of continuity of particle 
displacement. This yields 


dt dx 


(42) 


where £(x,0,f) = C(x)e ,(n,/ ’ me) is the normal component of wall displacement, directed into 
the wall, evaluated at the wall surface. It is assumed that all lined surfaces are nearly 
parallel to the duct axis of symmetry so that there are no high flow accelerations 
(particularly accelerations normal to the wall) in the lining region and so that the 
description of the lining displacement in terms of x is equivalent to a description in terms 
of the arc length along the wall. These assumptions are consistent with reality, and greatly 


simplify the lining model. Since vn - ~ it follows that with harmonic time dependence 


K d 

V = (1 -i — r -~ )»•« 

" n, ax 


(43) 


The relation between acoustic particle velocity and pressure is 


m a 

Tl r dx 


(44) 


The relation between pressure and acoustic velocity potential is provided by the acoustic 
Bernoulli equation of equation (9). Equation (44) can be rewritten 
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The boundary integral becomes 


-// Pr^v.ds = in ,/ f WAfods * / / 

5 S S ^ 


'ffvPrK^eWS - ±ffw Pr M± 


V 


dx 


A Pr M r 


dx 


dS 


( 46 ) 


The first two integrals on a boundary where acoustic treatment is present are easy to 
implement in the finite element formulation because only continuity of acoustic potential 
is required. The admittance, A, is assumed piecewise continuous. The third and fourth 
integrals have continuity problems because of A and d$/dx and are not compatible with 
the weak formulation. However, integration by parts can be performed to reduce the 
continuity requirement. This process begins with the observation that 


jjw Pr M,-hA f ,t]dS - 


dx 


dx [ 




Ap r M r 


d± 

dx 


(47) 


dS 


and 


Vs 31 




i rr 3 


dS = —ff 

V/ ^ 

'iff fr'rK] 


WA 

dx 


dS 




(48) 


dS 


Stokes’ Theorem for non-planar surfaces is 



where n is the normal to S and dl is the incremental line element on the curve C bounding 
S. The first integrals in equations (47) and (48) become 




±ff± 
n/j a* 


waM®. 

dx 


dS = — <£ 

V 


WAp 2 r M r r 2 -^ 

dx 


dl 


(50) 


( 51 ) 


In arriving at equations (50) and (51) advantage was taken of the fact that the surfaces on 
which acoustic lining is present are very nearly cylindrical. The bounding curve C is 
considered to consist of segments parallel to the duct axis and circular arcs bounding the 
lined region. The integrals on the line segments parallel to the duct axis cancel because of 
continuity considerations. The bounding circular arcs are chosen to be located just outside 
the lined region so that the admittance vanishes. With these arguments the line integrals 
discarded. The weak formulation for the boundary condition' on the acoustically treated 
surface can now be written 


-/ / p - in,/ / wa p^js 

s s 


ffWAplM&dS 


-nil"*'*' 

s 



(52) 


Equation (52) is in a form which is appropriate for application of standard finite element 
techniques to generate "boundary matrices" which are appended to the element stiffness 
matrices of elements whose outer boundaries represent acoustically treated surfaces. 



SOLUTION METHOD 


When the matching procedure is carried out and assembly is accomplished by 
standard finite element methods, the acoustic field is discretized in terms of nodal values 
of pressure in the pressure elements, nodal values of velocity potential in the potential 
elements, and amplitudes of the duct acoustic modes at the interfaces. The force integral l f 
is produced by the evaluation of the weighted residual volume integral 

l, " * r /// ™-fdV (53) 

V 


formed in the pressure elements in which the blade force distribution tf b } is defined at the 
element nodes [1]. 

The nodal values of pressure or velocity potential can be recovered from the modal 
amplitudes at the interfaces by post processing using the transformations described. The 
modal amplitudes are themselves useful information as they can be used to quantify the 
acoustic modes generated by the source and reflected from the fan inlet or and exhaust 
exits. All other aspects of the implementation of the finite element method are the same as 
described in connection with previous work [1,2, 3, 9]. 

The set of algebraic equations which arises from the finite element formulation is 
solved by using the frontal solution method of Irons [10], modified to deal with unsymmetric 
problems. Solutions have been obtained with over 20,000 degrees of freedom (nodes) with 
good success. 

GEOMETRY AND FLOW FIELD CALCULATIONS 

In Reference [1] it was assumed that the mean flow field is everywhere uniform. In 
addition to simplification of the acoustic field equations, this eliminated the complication 
of producing input data for a nonuniform mean flow field. In the present formulation the 
nonuniform flow field is required to be known. In previous studies of turbofan inlet noise 
radiation [2-4] it was noted that this data was obtained from a finite element formulation 
for incompressible potential flow based on the same mesh as the acoustic calculations. This 
potential flow code has been extended to the geometry of the simultaneous forward and aft 
radiation problem which now arises. The issue of the shear layer which is present at the 
boundary of the exhaust jet and the surrounding flow is not addressed in this study. Also not 
included in the model are energy and momentum jumps in the mean flow field which occur 
across the rotor. In the last case, a much more detailed analysis of the noise generation 
mechanism would be required to include this in a rigorous way. 

An automated mesh generation procedure which was developed for the inlet 
radiation has been extended to cover the forward and aft radiation geometry required in the 
new formulation. The mesh is constructed from input data which describes the nacelle and 
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centerbody. 

VERIFICATION OF THE COUPLING SCHEME 

The technique used to couple the pressure region and the velocity potential region 
requires verification. In order to do this, a uniform duct example has been considered in 
which the duct terminations have been made reflection free by use of a mode matching 
procedure not unlike the scheme described here for the coupling of the two regions [6]. 
Figure 5 shows the example with a region in the center of the duct in which the propagation 
is described by the pressure formulation. The remainder of the duct is represented by a 
velocity potential formulation. At the downstream end of the duct (x=0) a single acoustic 
mode propagating in the positive x direction is introduced. The reflected mode amplitudes 
at x=0 and the transmitted mode amplitudes at x=l are computed. The expected result is 
the absence of reflected modes and the presence at x=l of only the incident wave. This tests 
not only the scattering introduced by the coupling procedure but also that caused by the 
reflection free termination formulation. 

Table 1 provides results which are typical. In this case an angular mode m= 10 with 
the first radial mode n=l at nondimensional frequency x\ r = 20 is introduced at x=0 with 
unit pressure amplitude. There are three radial modes which are propagating. The table 
shows the first five incident, reflected and transmitted velocity potential modal coefficients. 
The largest reflected mode coefficient amplitude, which is the reflection of the incident 
mode, is .009 % of the incident mode coefficient amplitude. The transmitted modal 
coefficient corresponding to the introduced mode is very close to the result which would be 

ik* l 

obtained analytically by accounting for the phase shift introduced by the phase term e x ‘ . 
is the axial wave number for the first radial mode and 1 is the duct section length. 

A second example in which the third radial mode is introduced produces even better 
results. This is not unexpected because the effective wave length is longer than for the lower 
order modes and the discretization errors associated with mesh refinement should be less 
critical. The maximum reflected modal coefficient amplitude is .0006 % of the input 
amplitude and the complex modal coefficient for the transmitted mode is again very close 
to the analytical result. These calculations are shown in Table 2. 

The calculations described were based on a finite element mesh which was chosen 
to be refined enough to cope with the specified frequency. Numerical experiments have 
shown that five elements per wave length (accounting for Doppler effects) is reasonable and 
the results here are consistent with this rule of thumb. 

The conclusion to be drawn here is that the mechanics of the coupling procedure as 
described are sound, producing spurious scattering only to a level which can be related to 
discretization errors and errors introduced by using a coupling procedure and an anechoic 
termination model based on finite acoustic mode expansions. 
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Table 1 Verification of Coupling - First Radial Mode 
Angular Mode m = 10 
Frequency ti r = 20 
Flow Mach M = -0.30 
L/R = 0.4 


Mode 

Incident 

Coefficient 

Reflected 

Amplitude 

Transmitted 

Amplitude 

Transmitted 

Coefficient 

Transmitted 
Coefficient (Th.) 

1 

0.00 + i 0365(4) 

030(-5) 

0365(4) 

-0.168(-1) - i 0324(-l) 

-0.171(4) - i 0322(-l) 

2 


0.98(-8) 

0.42(-6) 



3 


0.47(-8) 

0.83 (-7) 




Table 2 Verification of Coupling - Third Radial Mode 
Angular Mode m = 10 
Frequency r| r = 20 
Flow Mach M = -0.30 
L/R = 0.4 


Mode 

Incident 

Coefficient 

Reflected 

Amplitude 

Transmitted 

Amplitude 

Transmitted 

Coefficient 

Transmitted 
Coefficient (Th.) 

1 


0.86(-8) 

0.19(-7) 



2 


0.21(-7) 

0.43(-6) 



3 

0.00 + i 0.422(4) 

0.60(-9) 

0.422(4) 

-0.409(-1) + i 0.101(-1) 

-0.409(-l) + i 0.104(-1) 

4 


0.26(-6) 

0.27(-8) 



5 


0.28(-7) 

0.43(-9) 
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ACOUSTIC LINING PARAMETERS 


In carrying out calculations to demonstrate the effect of acoustic treatment on the 
far field acoustic radiation it is necessary to choose suitable admittance values. This has 
been done by using a lining optimization code which uses a Simplex scheme to determine 
the optimum impedance for a uniform duct. Optimization can be done by seeking the 
maximum attenuation in a given mode (the equivalent to the "Cremer" optimum [11] when 
mean flow is present), or on the basis of power attenuation with several choices of modal 
amplitude distribution. The choices in the code are (a) The "Cremer" optimum; (b) specified 
modal amplitudes and optimization on the basis of transmitted acoustic power; (c) equal 
modal amplitudes and optimization on the basis of transmitted acoustic power; and (d) 
equal modal power and optimization based on transmitted acoustic power. In the present 
study we have used the "Cremer" optimum based on the first radial mode, two specified 
modal coefficients defined by the source model, and the equal modal amplitude and equal 
modal power assumptions for sample calculations. The duct has been taken as annular with 
the radius ratio at the fan, and the acoustic treatment on the inner and outer walls is the 
same. For the cases involving transmitted power the lining length has been taken as 60 % 
of the fan radius. The lining design is constrained by the requirement that the real part of 
the normalized admittance be greater than zero and less than one. This restriction on the 
maximum value of the admittance, which has been used to demonstrate the optimization 
capability, limits the effectiveness of the exhaust duct lining. The Mach number is taken as 
that at the fan face. The frequency and angular mode number are chosen according to the 
characteristics of the source. Table 3 provides the admittances for the inlet and exhaust 
ducts for the four cases. 


Table 3 Optimum Admittances' 
Inlet Mach = 0.4 
Exhaust Mach = 0.4 
Radius Ratio = 0.3 
Angular Mode m = 1 
Frequency -n r = 6.40 
Frequency r\ f = 6.44 
Propagating modes: 2 inlet and exhaust 
Constraint: 0.0 < Real[admittance] < 1.0 



Radial Mode 1 

Specified 

Equal 

Equal 



Coefficients 

Amplitudes 

Power 

Inlet 

0.22 + i 0.45 

0.68 + i 030 

0.10 + i 038 

0.18 + i 039 

Exhaust 

1.00 + i 0.82 

0.97 + i 035 

1.00 + i 0.98 

1.00 + i 0.22 
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COMPUTATIONAL RESULTS 


Example calculations will be shown for EGV interaction noise for a ducted fan. The 
principal feature which will be addressed is the attenuation predicted in the far field due 
to the insertion of acoustic treatment in the inlet and exhaust ducts of an inlet of simplified 
geometry. The geometry is for a model scale ducted fan with radius 0.311 m (1.02 ft). The 
blade chord is uniform at 0.052 m (0.17 ft). A nondimensional rotor angular velocity oft| r 
= 0.8 is for a case of subsonic tip speed. At a speed of sound of 344 m/sec (1128 ft/sec) this 
corresponds to a rotor rotational speed of 8448 RPM. The geometry of the simple nacelle 
is shown in Figure 6. The forward flight speed of the nacelle is M=0.3 and the flow entering 
and leaving the fan is M = 0.4. In the EGV case considered here there are eight blades on 
the fan and seven exit guide vanes. The interaction mode number is m=l. The 
nondimensional frequency at the fan plane is x\ f = 6.44 which produces two propagating 
modes in both the inlet and exhaust ducts. 

Calculations have been made for the case of untreated walls in the inlet and exhaust 
duct and for four cases of optimal acoustic treatment on the duct and centerbody in the inlet 
and exhaust ducts. As noted previously, the four methods of optimization depend on how 
the modal amplitudes are chosen. The classical, or "Cremer Optimum", is obtained by 
seeking the impedance which will maximize the attenuation in a particular mode. This 
occurs when two modal attenuations coalesce. In the present study the first radial mode is 
chosen and the optimization process maximizes the attenuation in the least attenuated 
mode. The remaining three optimal treatments depend on the choice of the mix of modal 
amplitudes and a maximum decrease of acoustic power in a specified duct length is the 
objective of the optimization. As a consequence of the matching of pressure and potential 
regions of the solution, the source acoustic modal amplitudes are available. These are used 
to determine the optimal impedance for the case of known amplitudes. The other two cases 
of equal modal amplitudes and equal modal acoustic power circumvent the necessity of 
knowledge of the modal structure of the source. The lining length in both the inlet and 
exhaust ducts is 60 % of the fan radius on the nacelle and 87 % of the fan radius on the 
centerbody. The calculations in this example will shed some light on the most appropriate 
way to choose the optimal acoustic treatment for this configuration. 

The mesh for the example calculation is shown in Figure 7. The most important 
feature to note is the outer boundary of the computational domain which is a circle with 
center offset in the direction of the flow approaching the inlet. This is a constant phase 
surface for an apparent acoustic source at the origin, consistent with the use of wave 
envelope elements in the outer reaches of the domain to economically introduce the 
reflection free condition at "infinity". 

Figure 8 is a plot showing contours of constant sound pressure level around the 
nacelle in a plane through the axis of symmetry for the case when there is no acoustic 
treatment. It shows the presence of two lobes in both the forward and aft radiated fields. 
In the region of 90 degrees to the axis of symmetry other lobes appear which are the result 
of diffraction and interference between the forward and aft radiated fields. The peak of the 
forward radiated noise is about 2 dB less than the aft radiated noise (this observation is 
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more clearly seen in Figure 10). The contours near the inlet tend to be slightly ragged due 
to postprocessing velocity potential data at the element nodes to obtain pressure. The 
pressure requires derivatives of the velocity potential and it is known that the weak 
formulation does not ensure continuity of derivatives on element boundaries. A superior 
postprocessing procedure computes the pressure at Gauss points within the elements where 
continuity is assured. This has not been done in this investigation because the contour 
plotting software requires nodal values of pressure. In Reference [9] the more accurate 
scheme was used and shown to generate better contours. It is interesting to note that the 
contours are more smooth in the wave envelope region where interpolation includes the 
wave structure. 

When acoustic treatment is present the radiated sound field is considerably altered. 
Figure 9 shows the result when the "Cremer" optimum is used. The optimization process 
indicates that the exhaust duct lining is only about 25 % as effective as the inlet lining (in 
terms of dB per unit length attenuation). While the contours in Figure 9 are not identified 
quantitatively, it is found that the acoustic field in the forward arc now has a distinct single 
lobe, and this lobe is down more than 10 dB relative to the aft radiated noise which has had 
its peak level reduced by about 7 dB due to acoustic treatment. The aft radiated noise now 
has a broad single peak. Quantitative observations are more easily deduced from Figure 10. 

Figure 10 is the polar directivity of the radiated acoustic field on a circle with radius 
of 10 fan radii. The five cases are shown on this plot with the largest Sound Pressure Level 
normalized to 100 dB. The scale level shown of 105.97 dB indicates the highest SPL, so that 
the actual level on any of the curves is obtained by adding 5.97 dB. The radiation pattern 
in the unlined case can be compared to Figure 8 to identify lobes and to quantify the levels. 
The characteristics of the level curves in Figure 9 can also be identified on Figure 10. Due 
to the low angular mode number, m=l , the acoustic field has its maximum levels at 
relatively low angles from the axis of symmetry. 

The four cases of duct acoustic treatment show generally similar effects on the 
acoustic far field. Both attenuation and shift in location of the lobes is a characteristic of all 
four cases. One could view attenuation as the decrease in SPL in easily identifiable lobes, 
regardless of angular shifts in the lobes. This would provide the most optimistic definition 
because it would obscure the fact that the new lobes might move into previously low SPL 
regions. The definition of attenuation adopted here will be the decrease of SPL at a specific 
polar angle. This is the most conservative definition because it is adversely affected by low 
SPL levels in the field for the untreated nacelle. For this fan configuration two choices of 
optimization philosophy appear to produce superior results for attenuation in the far field. 
The "Cremer" optimum and the equal modal amplitude scheme give generally good results 
in the entire field, with the attenuation being 10 dB or more for forward radiated noise and 
sideline attenuation in the same range except for isolated angular locations. In the case of 
aft radiated noise the attenuation is 7 or 8 dB except in the neighborhood of 150 degrees 
where an interference dip in the unattenuated field is present. The equal power choice does 
well in the forward region but does not do as well in the aft region. It is surprising that the 
specification of modal amplitudes in the optimization scheme is the least successful of the 
four choices. One can rationalize this by noting that the design procedure focuses on a 
section of an infinitely long lining. Scattering due to the leading and trailing edges of the 
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lining is not present and input modal amplitudes are not altered by this scattering. These 
results can not be taken as general because of the constrained range of the real part of the 
admittance, and because of the modal structure of the source, which has only two 
propagating radial modes, both of which are strongly affected by the "Cremer" procedure. 
In fact, the two modes tend to coalesce, accounting for the loss of the distinctness in the far 
field lobes in the directivity. It might be expected that if many modes are propagating 
different conclusions would be made. 

CONCLUSIONS 

The finite element modeling of ducted fan acoustic radiation has been extended to 
include a source which represents the loading on rotating or stationary blades. The loading 
is introduced by volumetric forces in the acoustic field equations written in a pressure 
formulation for a small region in which the flow is uniform. Away from this source region 
an acoustic velocity potential formulation is used which includes the effects of a nonuniform 
mean flow field due to inlet flow and forward flight effects. The two regions are coupled 
by using the transition between acoustic velocity potential and pressure which is available 
in an acoustic Bernoulli equation. At the interface between the regions eigenfunction 
expansions are used to express both pressure and acoustic velocity potential in terms of 
modal amplitude coefficients which become unknowns in the finite element formulation. 
The matching at the interface also requires the appending of the natural boundary 
conditions which are appropriate for the two regions. A consequence of the matching 
procedure is that the modal amplitudes are calculated as part of the solution, thus providing 
information on the duct modes generated by the source. IQiowledge of the modal amplitudes 
is useful for the design of optimum acoustic treatment. The mechanics of the matching 
procedure has been tested in an example calculation and it is found that scattering at the 
interface is well within other sources of computational error associated with the FEM. 

An additional extension of the FEM model is the provision for the acoustic treatment 
of the surfaces of the nacelle and centerbody. Designated acoustic admittances can be used 
in a model of a locally reacting acoustic lining in both the inlet duct and exhaust duct. 

Example calculations have been made for the case of EGV interaction noise to 
investigate the effectiveness of acoustic treatment which is specified on the basis of several 
forms of optimization criteria. It is found that an optimum admittance based on the 
"Cremer" optimization or by the choice of equal modal amplitudes provides the best far field 
performance. Because of the complex nature of the propagation and radiation, the result 
must be considered as specific to the configuration in this example. Calculation of the 
radiated field provides a useful supplement to acoustic design techniques based on infinite 
duct theory. 
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Figure 1. Idealized view of rotor/EG V imbedded in a nacelle with a 
centerbody . 
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Figure 2. View of the inlet and exhaust ducts showing the pressure 
region near the blade row and the potential region away 
from the blade row. 



Figure 3. Boundaries of the computational domain showing the center 
body, nacelle with acoustically treated areas, and far 
field boundary. 
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Figure 4. Interface of a pressure region and a potential region 

emphasizing the outward normals to the surfaces on which 
natural boundary conditions apply. 
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Figure 5. Annular computational region used to test coupling 
procedure for pressure and potential regions. 



Figure 6. Geometry of simple nacelle shape used in example 
calculations. 
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Figure 7. Mesh generated for example calculations 
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Figure 8. Contours of constant Sound Pressure Level around the 
nacelle in the case of EGV interaction noise with 8 
blades and 7 exit guide vanes creating the angular mode 
m=l. The blade tip Mach number is M=0.8. The inlet and 
exhaust ducts are unlined. 
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Figure 9. Contours of constant Sound Pressure Level around the 
nacelle in the case of EGV interaction noise with 8 
blades and 7 exit guide vanes creating the angular mode 
m=l. The blade tip Mach number is M=0.8. The inlet and 
exhaust ducts are lined with a "Cremer" optimum lining. 
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Figure 10. Polar directivity at 10 duct radii for EGV interaction 
noise. 8 blades and 7 exit guide vanes. Blade tip Mach 
number is 0.8. Amplitude normalized to 100 dB. Computed 
level is obtained by noting the scale level of 106 dB. 
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ABSTRACT 

A finite element model is created for the generation, . 
propagation, and radiation of steady, rotor alone noise and 
rotor and exit guide vane interaction noise of a ducted fan. 
In the case of rotor alone noise the acoustic source is 
represented by a rotating lifting line of thrust and torque 
dipoles distributed radially on the blade. For a specified 
number of blades, angular mode harmonic, and rotor 
angular velocity, the acoustic field is described in a 
cylindrical coor&nate system reduced to only the axial and 
radial directions. The blade loading is imposed on a radial 
line at the axial location of the fan. The blade loading is 
assumed to vary linearly from the hub to the tip. The blade 
tip loading is determined by a specified thrust requirement 
and the inflow velocity. In the case of interaction noise the 
acoustic source is a stationary lifting line of torque and 
thrust dipoles which represents the fluctuating lift on the 
exit guide vane created by the velocity deficit associated 
with wakes in the steady velocity field behind the rotor. The 
fan and exit guide vanes are imbedded in a duct which can 
be contoured to represent a realistic installation. In the 
configurations considered in the present study, emphasis is 
on ducted fans or ducted propellers for which the by-pass 
ratio is very large. In this case the usual assumption is made 
that the fan, or propeller, is operating in a mean flow 
environment which is uniform and the same as the forward 
flight velocity. The flow acceleration in the inlet, 
acceleration in the fan duct, and jet free shear layer are not 
accounted for in the present model. The model accounts for 
the noise generation process, the propagation through the 
inlet and fan duct, and the radiation to the near and far 
field. The major issue addressed in the computational 
examples is the relationship between the far field radiated 
Sound Pressure Level (SPL) and directivity and fan tip 
speed. 

INTRODUCTION 

Ultra high by-pass ratio turbo-fan engines and ducted 
or shrouded propellers are attractive from the standpoint of 
propulsive efficiency. In addition there are possible 
advantages to be gained in radiated noise levels due to the 
imbedding of the propeller or fan acoustic source within the 
nacelle or shroud. An unducted propeller generates an 
acoustic field which tends to produce high levels on the 
sideline, and therefore may create unacceptable noise 
levels in the interior of the aircraft. A ducted propeller is 
restricted in the way in which it can radiate to the near and 
far field. It is known that steady, rotor alone noise, created 
by blade loading, is a principle source mechanism for 
unducted propellers. It is generally assumed that in the case 
of a ducted propeller the rotor alone noise is not 
propagated to the far field if the tip speed does not exceed 


’the speed of sound. This result, due to pioneering work of 
Tyler and Sofrin (1], is true for rotor generated noise in a 
thin annulus with the absence of duct mean flow. 

In the case of ducted fans and propellers an additional 
source mechanism exists associated with the presence of exit 
guide vanes. The EGV operate in a helical velocity field 
behind the rotor, which is for the most part steady and 
defined in direction by the thrust of the rotor. The mainly 
steady character of the rotor generated velocity is 
periodically interrupted by the viscous wakes downstream of 
the individual rotating blades. The EGV produce lift in 
response to the rotor velocity field, and because of the 
fluctuating velocity field behind the rotor, produce 
fluctuating lift and provide an acoustic source mechanism. 
This interaction source mechanism was also addressed by 
Tyler and Sofrin [1] and was shown to have the potential for 
the creation of acoustic modes of very low angular order. 
Modes of this type will propagate and radiate with a 
directivity pattern which may produce high levels near the 
axis of symmetry. 

The purpose of the work presented here is to 
investigate the differences in the radiated acoustic fields of 
ducted and unducted propellers of the same thrust operating 
under similar conditions. Hanson [2] has created a 
comprehensive acoustic model for unducted propellers 
which accounts for spanwise andchordwise details of the 
blade loading. It is not the intent in the present study to 
focus on such a refined model. Instead, the approach is to 
generate a very simple source model, similar to the classic 
lifting line theory suggested by Gutin [3], to concentrate on 
the propagation and radiation effects introduced by the duct, 
and to compare the acoustic performance of similar ducted 
and unducted propellers based on the same source model. 

The finite element method (FEM) has been used in 
previous studies to model the wind tunnel acoustic testing of 
propellers and the free field acoustic radiation of propellers 
[4-7]. In the present study the FEM is used to model the 
ducted propeller in the free field. This combines the 
propeller modeling previously reported and some aspects of 
earlier work on the prediction of the radiated acoustic field 
from turbofan engine inlets [8-9]. 

The generation, propagation, and radiation of sound 
from a ducted fan is described in this study by the convected 
wave equation with volumetric body forces. Body forces are 
used to introduce the blade loading for rotating blades and 
stationary exit guide vanes (EGV). For an axisymmetric 
nacelle or shroud, the problem is formulated in cylindrical 
coordinates. For a specified angular harmonic the angular 
coordinate is eliminated and a two dimensional 
representation results. A finite element discretization based 
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Figure 1. 

Geometry of Shrouded Propeller or Rotor. 


where p* is the acoustic pressure, p 0 is the ambient density, 
c is the ambient speed of sound, and J* represents the body 
force per unit mass acting on the fluid. p Q J' is the body 
force per unit volume. Equation (1) is in dimensional form. 
In the development which follows a nondimensional form of 
equation (1) is used with the following scaling 


X = 


L 9 




t" is the dimensional time and x* is any of the linear spatial 
coordinates. The reference length L is the propeller radius 
R. The nondimensional form of the acoustic field equation 
is 

v. 

( 2 ) 


on nine node quadratic isoparametric elements is used. The* 
nacelle and center body or core engine are defined as rigid 
surfaces. The assumption is made that the bypass ratio is 
large enough so that the entire flow field is uniform, 
consistent with the usual model in propeller acoustic 
analyses. Features not modeled are the nonuniform flow in 
the inlet and fan duct, and the free shear layer in the fan 
duct exhaust jet. 

GgQm ttrY and Coordinate System 

In this investigation the acoustic field is conveniently 
represented in a cylindrical geometry with the axis of the 
propeller or rotor/nacelle designated as the x axis. It is 
assumed that the nacelle/centerbody combination is axially 
symmetric and that the inlet flow field is axially symmetric. 
Specifically excluded by this restriction are drooped inlets 
and nacelles for which the inlet duct is not circular. The 
acoustic field is not required to be axially symmetric, and it 
would be unlikely that it is. The acoustic field is periodic in 
the angular coordinate of the cylindrical system. It is 
represented as the components of a Fourier Series in the 
angular coordinate 0. The acoustic field for each angular 
component, or "angular mode," is represented by a field 
equation in only the axial and radial components x, r of the 
cylindrical system. 

Figure 1 shows an idealized geometry for a 
rotor/nacelle arrangement. Noise sources related to the 
rotating blades and the interaction of the blades with 
stationary exit guide varies can be modeled. 

The steady velocity field in and around the nacelle is 
assumed to be uniform. For many applications, notably the 
case of the ultra high bypass fan or the ducted propeller this 
is probably satisfactory. For other cases it may be necessary 
to model the flow in and around the nacelle. It may also be 
required to consider the effects of the shear layer in the 
interface between the fan exhaust and the surrounding 
steady flow. 

Mathematical Model 

The acoustic field is described by the converted wave 
equation with body forces 


*P' 


i py 

c 1 Of** 


p.v*7* 


( 1 ) 


The body force per unit mass J is related to the force 
exerted on the fluid by the rotating blade or stationary vane. 

The major deficiency in this model is the assumption 
that the interior flow and external flow are uniform and at 
the flight Mach number. This is required because a pressure 
formulation has been chosen to introduce the acoustic 
source model for the rotor or EGV via equivalent body 
forces acting on the fluid. This is consistent with previous 
models of propellers [4-7]. In the pressure formulation it is 
required that the flow field be uniform in order that the 
acoustic field equations can be reduced to the converted 
wave equation. 

By modeling the acoustic field using the pressure 
formulation some liberties have been taken. Even though 
for ducted propellers or high by-pass ratio ducted fan types 
of flow fields the internal and external flows over most of 
the region may be approximated as being uniform and 
tangent to the nacelle, a region will exist near the inlet lip 
where the flow is clearly not uniform. The formulation of 
the problem using the convected wave equation does not 
account for this. Although the primary convective effect of 
the mean flow is modeled, effects due to localized flow 
gradients will not be included. 

In addition to the requirement that the flow be 
uniform there is a more subtle restriction that is introduced 
by the natural boundary condition for the convected wave 
equation that would require 

(Vp - M 1 & J) • n - 0 (3) 

ebe 

if no forced boundary condition exists, /f is the outward 
normal on the nacelle surface. On the nacelle surface the 
rigid wall boundary condition requires that 

Vp • n » 0 

which is equivalent to specifying that the acoustic particle 
velocity normal to the surface vanishes. Equation (3) does 
not reduce to this condition except where f . /f =, o* This 
is clearly violated where the surface normal is not 
perpendicular to the duct axis, and particularly near the 
inlet lip. The natural boundary condition (3) is used here 
with the argument that its apparent failure where i ■ fi + Q 
is an artifact of the assumption of uniform flow. No 
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apparent effect of this approximation can be seen in the 
computational results. 

A second modeling option has been used in previous 
studies of acoustic radiation from turbofan inlets [8,9]. The 
assumption is made that the mean flow and acoustic field 
are irrotational. The acoustic source is introduced through 
the specification of acoustic modal amplitudes at the fan 
face. This is accomplished by including suitable boundary 
terms. It is therefore not required to include the body forces 
in the momentum equation. This allows a first integral of 
the linearized momentum equation to be obtained, which 
can be interpreted as an acoustic Bernoulli equation. The 
combination of the continuity equation and the acoustic 
Bernoulli equation written in terms of the acoustic velocity 
potential provides the basis of the mathematical model in 
[8,9]. 


The velocity potential formulation does not seem to 
be appropriate if the source model is to be introduced 
through equivalent body forces. The critical feature which is 
lost is the first integral of the momentum equation leading 
to the acoustic Bernoulli equatioa This integral is not 
possible unless the body force is derivable from a potential. 
This does not appear to be the case for the types of body 
forces required to represent a propeller, rotor, or EGV. 



Figure 2. Velocities and Inflow Angles 

Experienced by a Blade Element. 


The lift per unit span at the tip is 

, *“i pt ' ,[1 + (;)’ ] w*. (6) 


c(r) is the local blade chord and c t is the section lift 
coefficient taken as constant in the present study. The lift 
per unit span as a function of radius is 


The finite element formulation of the converted wave 
equation (2), is the same as described in previous work 
related to propeller acoustic radiation [4-7]. Nine node 
quadratic isoparametric elements are used with special 
attention given to the elements spanning the propeller and 
containing the source terms. A frontal solver is used and 
models with as many as 18000 degrees of freedom have 
been handled routinely. 

The models used for the propeller or rotor and exit 
guide vane acoustic sources are discussed in the following 
sections. 


Kr) - 4 



c(r) a(r) 


(7) 


For the present investigation it is assumed that the 
blade loading is linearly distributed from root to tip, that is 


Kr) = (8) 


Bl ade L oading 


The thrust on N B blades is given by 


The blade loading of the propeller or rotor will be 
considered as the only source of rotor alone noise. No 
effects of blade thickness will be modeled in this 
investigation. Blade loading will be based on isolated lifting 
surface theory using a strip analysis. The discussion of 
Dommasch, Sherby and Connolly [10] is directly relevant to 
the following development. 

Figure 2 shows an airfoil section at the radius r from 
the hub. The local angle of attack of the section at radius 
r depends on the inflow velocity, U, the relative velocity due 
to rotation rn/U, where n is the angular velocity of the 
rotor or propeller, and the blade twist <p . The velocity seen 
by the section is 
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The angle B defining the direction of the velocity seen by 
the section and the angle of attack a are given by 
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Rj is the inner radius of the propeller or rotor. The 
required tip loading for a given thrust can now be 
determined from equation (10) 


'a? 
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( 12 ) 
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The thrust and torque loading are 


I, = j* cosp (i) , l m ■ ^ stop (-£) (13) 

Other types of loading can be considered by 
reformulating equation (8) and the subsequent analysis. 
Since the investigation reported here centers primarily on 
comparisons of unshrouded and shrouded propellers, it is 
not deemed critical to precisely specify the loading. 

The lift distribution discussed here is dimensional, 
and thus based on full scale parameters. 

Fqiqt Alone Noise 

Rotor alone noise generation is viewed from the 
perspective of a volumetric force fixed in space which is 
active during the passage of a blade with its associated lift 
distribution. It is assumed that the duration of passage of 

the blade past a fixed point is t s — where a(r) is the 

projection of the blade chord on the rotor plane. The 
strength of the volumetric force representing the blade 
passage is taken as the negative of the lifting pressure 
differential across the blade, approximated . by l(r)/c(r) 
where c(r) is the local chord. 

The strength is 


p.7- "fo (y ‘ )? * (14) 

where y’ is a dimensional coordinate normal to the blade. 
The lift per unit span here is taken to act normal to the 
blade chord, which is oriented by the twist angle <p . The 
unit vector e M is taken normal to the blade. The 
relationship between x\ the dimensional axial coordinate, 
and y\ the normal to the blade chord, for 0 = constant, is 


y * a x’cos<|> 


By making use of the property of the Dirac delta function 
that <s(ax) = ^(x), equation (14) can be written 
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c(r)cos0 is the projection of the blade chord on the rotor 
plane so that 


p.7= - ac**)?. 06) 

At a fixed angular position 6 = 0, for the successive passage 
of N b blades 
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The body force per unit volume is periodic with period 


r * 


2k 


and can be expressed in the Fourier Series 

p.7» - 5(0 E c ^‘ v ( 18 ) 
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The body force per unit volume can be resolved into 
thrust and torque components on the basis of the angle 
which orients the velocity seen by the blade. The 
components are 

P./, = - & cosp 5(0 E (20) 

a »— 

P,/.»Msiape«i: (21) 

a 

Note that these body forces are on the fluid and are 
therefore opposite to the corresponding forces on the blade. 
The thrust component is rearward and the torque 
component is in the direction of blade rotation. 

The development to this point is for the reference 
location 0 = 0. At any other location the same temporal 
event occurs with the phase lag At = 0/n . Hence at any 
angle, in non-dimensional form 

— (-^) cosp 5(*) X c . eim>,,V 
pc* <* 

( 22 ) 

/. - -L (M) stop 5(x) i C m 
9*. a m~ 

(23) 

The body force per unit volume has frequencies which are 
integer multiples of the blade passgage frequency N B n. The 
nondimensional form is obtained by noting that the 
nondimensional body force per unit volume is obtained from 
the dimensional form by dividing by c 0 2 /R. Furthermore, 
the dimensional Dirac delta function with the dimensional 
argument is nondimensionalized by division by R. 


E GV I nteraction Noiss 

In order to estimate the noise generating mechanism 
of the exit guide vanes (EGV) in their interaction with the 
rotating blades, a simplified model has been constructed. 
This model assumes that the EGV are on the average under 
the influence of a steady lift dictated by the magnitude and 
direction of the absolute velocity field leaving the rotor. 
This can be approximated from a knowledge of the steady 
blade loading in the rotor stage. The flow field downstream 
of the rotor is not steady, but is interrupted by the wakes 
downstream of the blade trailing edges. The velocity deficit 
in the wake, which is dependent on the distance downstream 
of the blade trailing edge, creates a fluctuating lift on the 
EGV. It is this fluctuation which becomes the noise source. 
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r = RAOIUS OF BLAOE SECTION 
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span is 
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( 31 ) 


is the period of passage of the reference blade in front of 
the reference EGV. The duration of the reduced velocity 
is established by the wake thickness h according to 


Figure 3. Velocity Triangles for Flow 
Behind the Rotor. 


hT 

t = 

2itr 


( 32 ) 


The model used here is based on quasi-steady strip theory With this picture of the fluctuating lift, the Fourier 

aerodynamics and is intended only as an estimate of the series for this periodic event is 
actual source mechanism. 


From the elementary theory of axial flow 
turbomachines [11] the absolute inlet velocity, v a , the 
velocity relative to the rotor at the inlet, v the velocity 
relative to the rotor at the exit, v^, and the absolute velocity 
at the exit v 2 can be determined from the velocity triangle 
of Figure 3, The pressure rise across the rotor is given by 

Ap = j U 1 (tan 2 p, - tan J p } ) (24) 

The pressure rise is balanced by the thrust component of 
force per unit span 

N b l^r) - 2nrAp (25) 

Therefore, for a known thrust loading 

W = 2 U 1 (tan 2 p, - tan 2 Pj) (26) 

2 N b 

For a local blade section & x and fij are defined by 



where q* = 1/2 pU 2 is the approach dynamic pressure. 
From the velocity triangle 

a 2 = tan" 1 (-5l - tanPj) , v 2 = Usj\ + tan 2 a 2 (28) 

The steady lift per unit span on an EGV is estimated to be 

l , * \ pv 3 C, e (i «J (29) 

with the assumption that the EGV has no twist c e is the 
chord of the EGV and the lift curve slope is c^. The steady 
lift per unit span on the EGV can be estimated from 
equation (29) at the same time that the rotor blade loading 
is calculated. 


/„(') ■/,-(!- «)', £ F. «‘* a ' ( 33 ) 

»•— 

where 

sin nDt + i(cos nQr - 1) 


The fraction of the steady lift experienced by the 
EGV during the encounter with the blade wake is assumed 
to depend on the distance between the rotating blade and 
the EGV. The measure of this separation is the distance 
in blade chords of the EGV leading edge behind the 
rotating blade leading edge. The velocity downstream in the 
wake is given by 

■ (1 - 7 ) V, (35) 

x d 

The parameter e is therefore 



In equation (29), the steady lift is associated with the 
steady flow. In the acoustic formulation the fluctuating load 
is the acoustic source, so that body force per unit volume 
will be related to the fluctuating lift per unit span given by 

Kt) = -(1 - t)l, ± F.e^ (37) 

a— • 

The EGV is represented spatially by a lifting line on 
which acts the force per unit volume 

p/ - /(r) «(**) 8 (rQ) ( 38 ) 

The Dirac delta function is written as a periodic function 
with period 2rr, 

s m - » - E G .<" . G . = t- <? 9 ) 


The load fluctuations on a reference EGV are 
stipulated to occur due to the velocity deficit in the wakes 
behind the reference rotor blade. The unsteady lift per unit 
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The body force per unit volume, the reaction to the lifting 
force, is 

p/u - (1 - «)*, ^ £ F, c™ ± G.e" (40) 

r I— «• — 

Equation (3) represents the fluctuating body force per unit 
volume on the reference EGV, numbered vane I, due to the 
passage of the reference blade, designated blade 1. 

For other than the reference vane there are temporal 
and spatial phase shifts 


2 K 

k = 1. N, 

(41) 

1 ~ 
K ’ 

k » 1, N y 

(42) 


Equation (41) is the phase shift in time at the vane k due to 
the intervane angular spacing 2 jt/N v and the angular 
velocity of the blade n. Equation (42) is the spatial phase 
shift due to the angular spacing. Hence at vane k due to 
the passage of the reference blade 


Pfu -(l-Ol.-22-i £ £ F P.‘ 




At vane k there is a further temporal phase shift due 
to the passage of blade / 


(/ - l)Af s (f * 1) ““ (44) 

The body force per unit volume for this interaction is 

P 7 a . (1 - e)/ m £ £ f,G„ 

r a — (45) 

- (/-I) -£-) Me-tk-i) 

x * *•“ ^ 

The total acoustic source contribution is obtained by 
summing over all vanes and blades 


p/ = (l-e)/,-^£ £ £ £ F A G. e * a ‘ e > 


Equation (46) can be simplified by making use of the 
result that 


= 0 , n + tigN 

s N » » = n/f, n B » 0, *1, ±2, . . 
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In equation (46) for a fixed value of m and n, consider the 
summation over k and t. This summation will vanish unless 

(47) 

(45) 


m - njf v - n/f B (49) 

The case n B = 0 is not relevant acoustically because it 
contributes to the time invariant part of the body force per 
unit volume. With these observations the body force per 
unit volume can be written 


n m V'a » n B “ ±l > ±2 » • • • 
m + n » n//, n Y = 0, ±1, ±2, . . . 

Equation (48) can be replaced by 



sin/itf-Gr + i(cosn// B Qt - 1) i 

F = 2 2 1 , G => — 

■ 2nnN s * 2 k 

Since the EGV is assumed to have zero twist the 
angle of attack is given by a 2 . The body force per unit 
volume is then resolved into thrust and torque components 
according to the non-dimensional relations 


f, - ^ Rr sia “i a <W t £ F ' G ‘‘ M ' ,a ‘ 
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FINITE ELEMENT FORMULATION 


Equation (2) is the field equation which governs the 
radiated sound field generated by the distribution of body 
forces J which are defined in the case of rotor alone noise 
by equations (22) and (23) and in the case of EGV 
interaction noise by equations (51) and (52). References 
(4-9) give details of the finite element discretization of 
equation (2). Here the important features are reiterated. 

A Galerkin weighted residual formulation seeks a 
solution for the acoustic pressure p among the class of 
functions C l with continuous first derivatives which satisfy 
the weighted residual statement 


///>,[*• J r - f) - 2iTiwj£ * I ftdV 

-ft W, - M 1 & T) ■ n - <yp - M 2 & 0 • ;?Us 

J ax ax 

= 0 

(53) 



for piecewise continuous weighting functions Wj. Equation 
(53) is the combination of the volume weighted residual of 
the field equation and the surface area weighted residual of 
a combination of terms which turns out to be the natural 
boundary condition. The overbar on p signifies that it is to 
be designated on S. The weak form of equation (53) is 
obtained by using the divergence theorem to reduce the 
level of continuity required of the solution p. In this form 
a solution p is sought in the class of continuous functions C\ 
with piecewise continuous first derivatives such that 

VW, • - A/ 2 & <j + 2ir\MW t - 

■ // vw, •/ dV * If W t ^Vp - M* dS 

(») 



Figure 4 . 


Computational Domain for FEM Formulation. 



for all weighting functions Wj in C\ Use is made of the fact 
that J , the body force, vanishes on the boundary S of the 
domain V. The body forces are distributed over the volume 
V f . n is the unit normal vector out of the volume V. 

Because of the axial symmetry of the geometry of the 
nacelle, equation (54) is implemented in cylindrical 
coordinates. Figure 4 shows the bounding surfaces of the 
computational domain in a 0 = constant slice. The surfaces 
of the center body S c and the nacelle S N are rigid so that 
Vp • n » 0- On these surfaces the natural boundary 
condition is given by the surface integral in equation (54), 
which if it were to vanish on S c and S N would require that 

(Vp - M 1 & o • n = 0 (55) 

dx 

on S c and S N . As noted previously, the pressure formulation 
used here does not allow the specification of the details of 
the mean flow in the neighborhood of the nacelle lip and 
the centerbody tip due to the restriction of a uniform axially 
directed flow inside and outside the nacelle. The important 
features that are missing are the tangency of the mean flow 
in the vicinity of the nacelle lip and centerbody tip and the 
flow gradients which exist there. The flow gradients exist 
over a length scale which is short compared to* a wave 
length and can reasonably be neglected. The boundary 
condition of equation (55) is derived on the basis of the 
restriction of uniform axial flow and hence that the normal ff 
is perpendicular to the x axis so that f * n - 0- Just as the 
flow is not tangent to the surface of the nacelle near the 
nacelle lip and centerbody tip, T * n does not vanish. In the 
pressure formulation described here, equation (55) is taken 
as the natural boundary condition for all rigid surfaces with 
the argument that in a more exact formulation, for example, 
the velocity potential formulation of References (8,9], the 
equivalent boundary condition would be based on the 
tangential component of mean flow, and would therefore 
vanish. The validity of these assumptions can only be tested 
by checking the computational results and looking for 
anomalies in the radiated acoustic field near the nacelle lip 
and centerbody tip. 

In the case when the propeller is unshrouded, the 
surface integral on S, and does not exist (in the present 
study the center body does not exist if the propeller is 
unshrouded). For the unshrouded case the uniform axial 
flow assumption is also not rigorously true because of flow 
gradients which must exist to account for the momentum 
increase which creates the thrust No known propeller 


acoustic model accounts for this. In this context, the 
formulation for the shrouded case is viewed as equivalent. 
With these arguments, equation (54) is taken as the weak 
form with the surface integral existing only on S*, the 
reflection free outer boundary. The important features of 
the finite element mesh used in the discretization of 
equation (54) are discussed in [4,8,9]. 

The volume integral over the distribution of body 
forces representing rotor alone or EGV interaction noise 
generation requires some special treatment because of the 
nature of the loading in the form of a Dirac delta function. 
Details are discussed in Reference [4]. 

The large set of algebraic equations created by the 
FEM discretization of equation (54) is solved by using a 
frontal solution routine. The resulting nodal pressures are 
postprocessed to create an acoustic directivity pattern, which 
is a contour of equal sound pressure level (SPL) on an x-r 
plane sliced out of the cartesian system, side line SPL, which 
is a plot of SPL on a line parallel to the axis and at a 
specified distance, and polar SPL, which is a plot of SPL at 
a fixed polar radius centered on the nacelle. Polar SPL 
results for several cases can be superposed on a summary 
polar SPL plot. Only polar SPL results are shown here. 

COMPUTATIONAL RESULTS 

In this study both rotor alone noise and EGV 
interaction noise for a shrouded propeller will be compared 
to unshrouded propeller noise for a fixed thrust. In the 
rotor alone case there is a distinct difference in the subsonic 
and supersonic tip speed cases, while in the EGV inter- 
action no fundamental difference is attributed to the 
difference in rotor speed. 

A Nacelle and Propeller Configuration 

A model scale propeller and nacelle is considered 
here. The propeller has four or eight blades and is of 
dimensional radius 0.311 m (1.02 ft). The blade chord is 
taken to be uniform at 0.052 m (0.17 ft). The non- 
dimensional propeller angular velocity is taken as q 3 0.8 
and tj = 0.9 in the subsonic case and ^ = 12 in the 
supersonic case. For a speed of sound of 1125 ft/sec; this 
corresponds to angular speeds fl « 8426 RPM and n = 9479 
RPM for the subsonic case and n » 12639 RPM for the 
supersonic case. The nacelle geometry is shown in Figure 
5. In the unshrouded propeller case no centerbody is 
present This has only a slight effect on the propeller 
loading. The source location, whether rotor alone or EGV 
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Figure 5. 

Idealized Shrouded Rotor for Example Cases. 


interaction, is just ahead of the center of the nacelle. The 
flow velocity inside and outside of the nacelle is M = 0.4. 

B. Rotor Alone Noise 

Figures 6 and 7 are summaries of the polar radiation 
directivity for four and eight blade shrouded propellers with 
comparisons to similar unshrouded propellers. Both 
supersonic and subsonic tip speeds are shown. Figure 6 is 
the four blade case. For supersonic tip speed, ri = 1.2, the 
(4, 1) mode (fourth angular, first radial) is cut off with cutoff 
ratio $ 4l (12)= 0.988 while at subsonic tip speeds it is cut 
off with cut off ratios f 4l (0.9) = 0.74 and £ 41 (0.8) = 0.66. 
The corresponding attenuations based on the cut off ratio in 
the duct length of 1.3R are 10 dB, 44 dB, and 49 dB. It is 
seen that the peak radiated noise at a distance of ten duct 
radii is 120 dB created by the unshrouded propeller (this 
sets the scale level for the plot). The shrouded supersonic 
propeller has a peak level of about 109. The unshrouded 
subsonic propeller with ^ = 0.9 peaks at 114 dB and the 
case with r\ - 0.8 peaks at 112 dB. The corresponding 
shrouded propellers peak at 80 dB and 74 dB. The 
relationship of the SPL levels of the shrouded propellers 
below those of the unshrouded propellers is consistent with 
the projected attenuations due to the cutoff phenomenon, 
though not numerically equivalent The comparison of the 
SPL levels of the two subsonic shrouded propellers with the 
supersonic shrouded propeller shows additional attenuations 
for the subsonic cases which are also consistent 

For the eight blade propeller a somewhat different 
picture emerges, as shown in Figure 7. For supersonic tip 
speed the (8,1) mode is propagating with l %l (1.2)= 1.085, 
but for the subsonic tip speeds it is cut off with £ 8l (0.9) = 

0.81 and (0.8) = 0.72 with calculated attenuations of 70 
dB and 83 dB. Reference to Figure 7 shows that the 
shrouded propeller at supersonic tip speed creates the 
highest SPL and sets the scale level at 130 dB. The 
corresponding unshrouded propeller has a peak level about 
10 dB lower. The unshrouded subsonic propeller at q = 
0.9 peaks at about 112 dB and the corresponding shrouded 
case is heavily attenuated at only about 74 dB. This 
attenuation of 56 dB with respect to the shrouded 
supersonic propeller is consistent with the cutoff 
phenomenon. A similar result applies when ^ = 0.8. The 
unshrouded propeller peaks at 108 dB while its shrouded 
counterpan peaks at 39 dB (hardly visible on the figure). 


This attenuation of 91 dB with respect to the shrouded 
supersonic propeller is also consistent with the cutoff 
phenomenon. The interesting feature here is the high peak 
SPL of the shrouded propeller. The (8,1) mode for the 
supersonic case is cut on, as opposed to the (4,1) mode 
being slightly cut off for the four blade propeller. It appears 
that the mechanics of wave propagation in the duct 
enhances the radiation of the ducted propeller noise when 
the mode is cut on. 

C. BGV Interaction Noise 

To demonstrate the radiation of EGV noise a case 
with eight rotating blades and seven stationary vanes located 
one blade chord downstream is considered. The lowest 
angular mode excited at blade passage frequency is M<, = 1. 
Figure 8 shows the radiated sound field for supersonic and 
subsonic rotors at r\ = 1-2, n 53 0.9, and r\ = 0.8. The 
radiation patterns are similar with high SPL near the axis of 
symmetry, characteristic of the well cut on M* = 1 mode. 
There is also significant SPL away from the axis due to the 
presence of higher order radial modes which are 
propagating. An interesting feature is the relative 
insensitivity of the peak SPL to the rotor speed, which is in 
contrast to the great sensitivity observed in the rotor alone 
noise. It is difficult to confidently compare the peak SPL in 
the rotor alone and interaction cases, since both source 
models are approximate and are probably not entirely 
consistent. It is safe to note that the peak SPL in the EGV 
case is much higher than in the case of subsonic rotor alone 
noise. EGV noise clearly becomes dominant in the subsonic 
case. 

SUMMARY 

Several important observations can be made. 

1) Contrary to the usual understanding of the Tyler and 
Sofrin result [1], supersonic tip speed rotor noise can be cut 
off if the tip Mach number is only slightly in excess of unity 
and if the number of blades is relatively small. If there are 
many blades, the fundamental angular mode number is 
large, and the Tyler and Sofrin result for thin annuli 
becomes more relevant. 2) Shrouding of subsonic tip speed 
propellers is a very effective means of controlling rotor 
alone noise. 3) There appears to be no benefit in terms of 
the peak radiated SPL for shrouded supersonic propellers 
when the fundamental mode is propagating. 4) For 
shrouded subsonic rotors, EGV noise becomes the dominant 
source. 
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APPENDIX B. MESH GENERATION SCHEME 

This section discusses the generation of the finite element mesh which plays a 
vital role in the formulation of the problem. Figure 1 shows an eight- node parent 
element with the local numbering of its nodes. For the convenience of constructing 
the mesh, the entire computational domain has been divided into three regions. 
Figure 2 illustrates the three regions clearly. Region I occupies the interior of 
the nacelle, both forward and aft of the propeller, region II extends from outside 
the inlet to the transition boundary C\ and region III (the wave envelope region) 
extends from C\ to the outer boundary (?«, • 

A. REGION I 

Due to the complex nature of the acoustic field inside the nacelle, a fine mesh 
is generated in order to resolve the variation in acoustic properties. It is separated 
from region II by a circle which we shall call the highlight circle. There are, in fact, 
two highlight circles, one fore and the other aft of the propeller. The highlight 
circles axe drawn from the nacelle tip (also known as the highlight) in such a way 
that their centers lie at the point of intersection of the x-axis and a line passing 
through the tip of the nacelle at 45° with the x-axis (see Figure 3). 

The inner surface of the nacelle C n extends from the propeller to the tip of 
the nacelle, both fore and aft. The centerline of the inlet geometry extends from 
the intersection of the centerbody curve with the x-axis (fore and aft) and the 
intersection of the highlight circles with the x-axis (fore and aft respectively). 
Three-node quadratic line elements lie along the inner surface of the nacelle, the 
centerbody and the centerline. The coordinates of these nodes are given as input 
to generate the mesh in region I. The number of input nodes on the inner surface 
of the nacelle and on the centerbody and centerline is the same, to produce a mesh 
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of generally rectangular elements. 

The vertical line segment x = 0 between the inner surface of the nacelle and 
the centerbody has also been divided into several elements not necessarily of equal 
width, each to be represented by a three-node quadratic line element. Figure 4 
illustrates the meshing scheme in this region. The “vertical” element boundaries 
inside the nacelle are formed by arcs of circles. These arcs are drawn through 
corresponding nodal points on the upper and lower boundaries (for example, the 
fifth node on the nacelle inner surface and on the centerbody and centerline, 
counted from the line x = 0) with the center of the circles lying on the x-axis. 
Such circles are easily constructed els illustrated by Figure 5. ( x\ , jq) and (x 2 , 
y 2 ) are coordinates of any two corresponding nodal points on the nacelle inner 
surface and the centerbody and centerline respectively. Then the x-coordinate of 
the center of a circle passing through these two points and having its center on 
the x-axis is given by 

(x| -x\ + yl-y\) 

*• — (1) 

Now, to preserve the rectangular mesh, each of these circular arcs should have 
the same number of three-node line elements on them and this should equal the 
number of three-node line elements on the line segment corresponding to x = 0. 
Therefore, each of these arcs is divided into the same number of elements with 
the same fractional length (fraction based on the arc length) as on the vertical r 
(x = 0) axis. Thus, the nodal coordinates of the line elements on each of these 
circular axes is determined. The vertical column of elements adjacent to the right 
of x = 0 is chosen to represent the propeller. The nodal coordinates are stored 
in a topology array AD(I, J, K) where I is the global element number, J is the 
local node number and K = 1 assigns the x coordinate, K = 2 the r coordinate, 
respectively, to the array. Proper connectivity relating the local node to its global 
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numbering is also generated and stored in a connectivity array AN(I,J), where 
I = element number, J — local node number. This array stores the global node 
number for each node. 

The global node numbering goes from top to bottom of each of the circular 
arcs starting from x = 0 and alternating between the fore and the aft of the 
propeller (see Figure 6). The element numbering is down each column of elements 
between adjacent circular arcs sequencing from x = 0 and alternating between 
fore and aft of the propeller (see figure 6). 

B. REGION II 

The mesh in region II becomes polar in nature essentially because of the 
configuration of the domain outside the inlet duct. Region II is separated from the 
wave envelope region III by a constant phase circle, as described previously, whose 
x-intercept is given as input. The constant phase circles are expanding radially 
with the local speed of sound (c) and their centers are moving away along the x- 
axis with the speed of uniform exterior flow (U) (see Figure 7). This phenomenon 
is very similar to the successive circles of outward ripples created on the surface of 
still water when a pebble is thrown into it. The only difference is that in still water 
the centers of the successive layers of outward moving circular ripples coincide and 
here the centers of the constant phase circles move at a constant velocity. 

From Figure 7, we obtain the equation of a constant phase circle (of radius 
R e ) displaced along the positive x-axis with velocity U (positive direction of U is 
indicated in Figure 7) 


(x - Ut) 2 + r 2 = R 2 


( 2 ) 
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where R e = ct is the radius of the circle Therefore, 

(x - U—Y + r 2 = R] (3) 

c 

or, 

(x — MRcY + r 2 = R* (4) 

where M is the Mach number of the uniform exterior flow. By setting r = 0 in 
equation (4) we obtain the x-intercept of the circle 

x = (M ± l)i2 e (5) 

The positive sign corresponds to the x-intercept on the positive x-axis, 

2 = (1 4- M)R e 

while the negative x-axis corresponds to the x-intercept on the negative x-axis, 

x = —(1 — M)R C 

The circles can be expressed in polar coordinates R and $ by, 

(. RcosO - MRcY + ( RsindY = R\ (6) 

Solving for R in terms of 0 yields, 

R = — M 2 + (Mcosfl) 2 - McosO] (7) 

Hence, the radial distance R at every angular position 6 on the outer boundary of 
region H is known. 

The outer surface of the nacelle which forms a part of the inner boundary of 
region II has three-node quadratic line elements along it. Since the mesh generation 
in region I precedes that in this region, the coordinates of the three-node line 
elements lying along the two highlight circle arcs are known. The nodal points 
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on the outer surface of the nacelle and those on the highlight circle arcs serve a s 
input for the mesh generation in region II (see Figure 8). 

In this region and also in the subsequent region III, the nodes have been 
generated on and along the acoustic rays from origin. Since the mesh is polar, 
the angular thickness of the elements increases with radial distance because the 
acoustic rays are radial lines diverging from the origin. To maintain proper aspect 
ratio of the elements in this region, the radial thickness of the elements should 
also increase accordingly along acoustic rays moving away from the origin. Now, 
corresponding to each nodal point on the outer nacelle surface and highlight circle 
arc, an acoustic ray is defined and its point of intersection with the outer bounding 
circle of region II is calculated (see Figure 8). Therefore, the radial distance along 
that ray in region II is known. This radial distance is then divided according to 
the number of elements required along the general direction of noise propagation, 
in geometric progression, from the inner boundary to the boundary C|. From 
elementary mathematics, we know that if ri, r 2 , . . . , r„ are n members of a series in 
geometric progression, then the members are related to each other in the following 
way, 

r 2 = cri ' 

rs = cr 2 

> 

r n = cr n _ x , 

where c is the common ratio of increment. So, the last member is related to the 
first member by 

r„ = c"- 1 ^ (9) 

Referring to Figure 9 where an acoustic ray intersects with the two boundaries 
of region II, it is obvious that the first and the last members of the geometric 
progression series, i.e. intersection points on the inner boundary and the outer 
bounding circle C\ respectively, are known. Since the number of elements n in the 
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radial direction of region II is an input, the common ratio of geometric progression 
is found out using equation (9) , 

( outer bounding circle \ " 

■: : : — — I 

inner bounding circle J 

Once the common ratio is known, the successive intervals are found out by multi- 
plication with the common ratio as in equations (8). Hence, the nodal points of the 
line elements along that acoustic ray are located. Geometric progression provides 
a gradual increment in the radial thicknesses of elements which is sufficient to a 
maintain proper aspect ratio. 

The nodal coordinate values are stored in rectangular cartesian form in a 
topology array AD(I, J, K) as mentioned before. The connectivity array AN(I , J) 
is also created. The element numbering and the global node numbering in this 
region is illustrated in Figure 6. 


C, REGION III 


Region III which consists entirely of wave envelope elements is bounded by 
the transition boundary C it the outer boundary Coo, and the x-axis. The wave 
envelope elements, as discussed before, are large elements bounded by acoustic 
rays and constant phase circles. The string of elements between any two successive 
constant phase circles is referred to as a wave envelope layer. The input for mesh 
generation in this region is the number of wave envelope layers and the x-intercept 
of the constant phase circles bounding each layer. Using equations (5) and (7), 
the inner and outer radii of the constant phase circles bounding each such layer 
is determined. The mean radius of each layer, which is just the average of the 
inner and outer radii, is also calculated. Since the mesh generation in region II 
is complete at this stage, the three-node line elements (note that a three-node 
line element forms a side of an eight or nine-node isoparametric element) on the 
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outer bounding circle C\ of region II have been located completely and their global 
numbering is also known. Therefore, corresponding to each nodal point on C\ t an 
acoustic ray is defined (see Figure 8) and thereby its points of intersection with 
the inner, mean and outer radii of each wave envelope layer are calculated. The 
rectangular cartesian coordinate values of these intersection points on which the 
nodes lie are stored in the topology array AD(I,J,K). The connectivity array 
AN (I, J) is similarly calculated as in region II. The element and the global node 
numbering follows after region II and is similar to region II. Since the mesh in 
region II is quite fine and that in region III is coarse, care should be taken to make 
a gradual transition in the size of the elements. 

After the mesh is generated, a connectivity check is performed to ensure a 
proper connection between local and global numbering of nodes and uniqueness 
of nodal coordinate values. 

E. SOME COMMENTS ABOUT THE FINITE ELEMENT MESH 

The acoustic radiation problem is highly mesh dependent but the mean flow 
problem is not very sensitive to the mesh parameters. Since both of these problems 
have been solved on the same finite element mesh, a mesh conforming to the 
acoustic parameters is desired. One of the important factors governing the mesh 
is the number of elements per wavelength which must always be maintained above 
a minimum value in the main direction of noise propagation to resolve the fine 
variation in acoustic properties. According to the rule of thumb the minimum 
ratio of the number of elements to the wavelength for quadratic elements should 
be 4 or 5. Here a somewhat crude estimate has been made to evaluate that ratio 
along the main direction of sound propagation. 

Since the nondimensional input frequency rj r (it is an input to the problem) 
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of the sound source on Cj is known, we obtain a ratio of the effective wavelength 
A e (the wavelength of the sound radiated is altered in the presence of mean flow) 
to the reference duct radius R in the following way : 

uR 2n R 

Vr = = — r— 

c A 

since 

A e = (1 — M)X 


therefore, 


Vr = 


2nR 

A e 


(1 — M) 


or, 


Ae 

R 


— (1 — M) 

Vr 


( 10 ) 


The Mach number M is positive if directed towards the inlet. Now if the number 
of elements per duct radius length is Nr and A is the average width of an element 
within that length, then 


R 

A 


= N r 


Therefore, using equation (10), the ratio of the number of elements per unit duct 
radius can be expressed as 


Nr 


Vr 

27r(l — M) 


Nx . 


( 11 ) 


where N\ t (= A e /A) is the number of elements per effective wavelength. For a 
specified number of elements per effective wavelength (for the elements used here 
N\ a is the goal), equation (12) can be used to determine the number of elements 
per unit of nondimensional length required. This varies as the flow towards the 
inlet varies, and would generally be highest within the nacelle near the propeller. 
The number of elements in the transverse direction within the nacelle or in the 
angular direction in region II is not as critical and is adjusted to maintain the 
aspect ratio of the elements. 
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Since the position of the constant phase circles bordering the wave envelope 
layers are user input, care should be maintained to make a gradual transition 
from the small conventional finite elements to the relatively large wave envelope 
elements. For this, the user should be aware of the radial thickness of the last layer 
of conventional elements along C\. Since the radial thicknesses of the elements in 
region II have been incremented in geometric progression, the radial thickness of 
the last element in region II on the x-axis is 

A r = (c n - c n ~ l )r 0 

where c is the common ratio of increment, n is the number of elements radially in 
region II and r 0 is the x-intercept of C\. This information is valuable to the user 
for making a smooth transition from region II to III. 
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Figure 2: Mesh generation regions 



Figure 3: Geometry of region I 





Figure 4: Node generation in region I 
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Figure 6: Element and node numbering in the finite element mesh 
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Figure 8: Node generation in regions II and III 




Figure 10: The computational domain 
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APPENDIX C. MEAN FLOW CALCULATIONS 

A. DERIVATION OF THE FLOW PROBLEMS 

Under the assumption of sufficiently low flow Mach number, the flow around 
the inlet of the shrouded propeller can be approximated as an incompressible one, 
and the equation for mean flow can be approximated by 

V 2 ^ = 0 (1) 

which is the Laplace equation. 

In Figure 10, the curves T in the x-r plane correspond to surfaces in the 
axisymmetric space around the inlet. The axis of symmetry r„ is not a physical 
boundary of the domain. In the axisymmetric integral formulation of the problem, 
the boundary integral corresponding to this axis (r = 0) vanishes. Therefore, no 
boundary condition needs to be specified. The far field boundary Too is circular 
(spherical in axisymmetric space). Though it is several duct radii from the inlet, 
the flow effects due to the presence of the inlet cannot be assumed negligible. 
The boundary condition on this curve will be discussed later. The nacelle r n 
and centerbody T c are impervious to flow. The curves T and T / 2 represent exit 
planes (convenient boundaries near the propeller where mean flow velocities are 
specified) fore and aft of the propeller, respectively. 

Since the differential equation (1) is linear, it can be split up into four different 
problems, each of which can be solved separately and upon employing the method 
of superposition, the velocity potential of the actual flow field can be determined. 
The velocity potential for the mean flow is decomposed as follows: 

4>o — <f>u + <j>p (2) 


where <f> u is the flow field due to the external uniform flow field only (without the 
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presence of the inlet) and fa is the flow perturbation due to the presence of the 
inlet only. The boundary condition to be applied at the boundarie Too is not clear 
until and unless we formulate the problem in terms of flow perturbation. Our aim 
is to formulate the entire problem in such a way that the fan face and the external 
flow velocity do not become dependent on the perturbations; rather they govern 
it. 


The perturbation velocity potential fa is further decomposed into 

4>p = <I>1 + fa + fa (3) 

where <f > i is the perturbation potential due to flow to a blank inlet (effect of the 
presence of the shroud in the external uniform flow), fa is the perturbation po- 
tential due to inlet flow through fore exit plane alone, and fa is the perturbation 
potential due to inlet flow through aft exit plane only. Therefore, the four flow 
problems may be posed as 


1. Problem I This problem represents the perturbation potential field due 
to a flow to a blank inlet. 


V 2 fa = 0 in n 

(4a) 

V<£i«n = — V<£ u • n on T/ 

(46) 

Vfa • n = —Vfa • n on T„ and r c 

(4c) 

_ . At 

Vfa • n = — r • n on 

M 


where fa is the external uniform flow velocity potential, Ai is a constant to be 
determined, n is the unit outward normal on the boundary and r is the outward 
radial vector on the outer boundary Too as shown in Figure 10. It is assumed that 
on the outer boundary Too the effect of the flow field is that of a simple source 
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placed at the origin. Hence, the velocity perturbation at the outer boundary is 
assumed to be radially directed inwards and inversely proportional to the square 
of the radial distance from the origin 1 . 


2. Problem II This problem represents the perturbation potential field due 
to inlet flow through fore exit plane alone. 


V 2 4> 2 =0 in fl 

(5a) 

V<f > 2 •n = Uf x on T f x 

(56) 

V<f> 2 • n = 0 on T„ and r e 

(5c) 

. . , -A 2 _ 

V<p 2 • n = r • n on Too 

R 1 

M 


where U/ x is the uniform exit face velocity, A 2 is a constant to be determined, n 
is the unit outward normal on the boundary and r is the outward radial vector 
on the outer boundary Too as shown in Figure 10. Here also the flow at the outer 
boundary is assumed to be that of a simple source placed at the origin. The flow 
perturbation is assumed to be radially outwards and varying as 1/R 2 , where R is 
the radial distance from the origin. 


1. Problem III This problem represents the perturbation potential field due 
to inlet flow through aft exit plane alone. 


V 2 $3 = 0 in fl 

(6a) 

V (f)^ • n = —Uf 3 on T/, 

(66) 

V<^ 3 • n = 0 on r„ and r c 

(6c) 


x The velocity field of a simple source varies as 1/R 2 . 
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V<^ 3 * n = — r«n on Too (6<f) 

where C7/ a is the uniform fan face velocity, Ai is a constant to be determined, n is 
the unit outward normal on the boundary and r is the outward radial vector on 
the outer boundary Too as shown in Figure 10. Here again the flow at the outer 
boundary is assumed to be that of a simple source placed at the origin. The flow 
perturbation is assumed to be radially outwards and varying as 1/R 2 , where R is 
the radial distance from the origin. 


3. Problem fV The uniform external flow field is generated by 


G 

cl 

• H 

O 

II 

3 

"O- 

> 

(7a) 

V<f> u • n = U 0 i • n on Too and Tj 

(76) 


where U 0 is the external uniform flow velocity. 

Problems I through III are boundary value problems with Neumann boundary 
conditions. Solutions to these problems are non-unique if the value of the unknown 
variable is not specified at one point in the domain fl. The problems also have 
' to satisfy a compatibility criterion which balances the flux of flow across different 
boundaries. This criterion fixes the values of the constants A\ through .A3 relative 
to the flow parameters and, therefore, they are not arbitrary. 

A linear superposition of the problems I through III gives us the overall bound- 
ary value problem of the mean flow 

V 2 <j> 0 = 0 in Q (8a) 

V(f> 0 • n = U/ on f/ 


V<f> 0 • n = 0 on T„ and r c 


( 86 ) 

(8c) 
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V<£„ • n = j±t • n - • n + ^fr • n 4- U a i • n on r M (8 d) 

The flow perturbation effects of the inlet at the outer boundary cure small 
due to the distance of the boundary from the inlet. Also it is to be noted 
that the perturbation boundary condition at the outer boundary Too for prob- 
lems I through III tend to balance each other. Therefore, under these conditions 
the superposed boundary condition (8d) on Too can be written approximately as 
V<£ 0 • n « U 0 i • n. The superposition of the elementary solutions from problems 
I, II, III, and IV is based on the assumption that the outer boundary condition is 
imposed at a large distance from the inlet. This effectively makes U / and U 0 inde- 
pendent of each other. For a given value of U 0 , any value of U/ can be chosen once 
the elementary solutions are available. Variations in U Q requires new potential 
flow solutions to be computed because the mesh depends on U Q . 


B. FINITE ELEMENT FORMULATION 


1. Flow to a blank inlet Since the mean flow field is axisymmetric, there is 
no variation of flow variables in the angular direction. Therefore, the test and 
trial functions are independent of the angular coordinate 6. Let rp be a real valued 
smooth function defined in the axisymmetric domain fi. Multiplication of the 
Laplace equation with the test function ip and integration over the domain yields 

[ VVilMn = 0 (9) 

Jet 

By using Green’s theorem, it is determined that 


f V<p\ • Vip dCl = [ rpV<pi»ndS 
Jn Js 


( 10 ) 


where S denotes surfaces of the axisymmetric volume fi. Since the problem is 
independent of 0, the volume integral becomes a surface integral in the x-r plane, 
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and the surface integral becomes a line integral, so that 


//< 


d<j>\ dxj) dfii dip 


) rdxdr = J 0V<£i • n rdT 


( 11 ) 


dx dx dr dr 

Incorporation of the natural boundary conditions into equation (62) results in the 
weak form of the problem 


//,< 


d<j > i dtl) d<f > i di\) 


) rdxdr — - J %l)V<f> 0 • n rdr + Ai f • n rdT (12) 

* Jroo it 


dx dx dr dr 

Therefore the following weak problem may be posed: find <j>\ : fi — > R 2 9 equa- 
tion (63) holds V smooth rj ) : H — ► R 2 . It is to be noted that <j>i and are suitable 
classes of functions whose derivatives are square integrable (from the space H 1 ) . 

A standard Galerkin finite element approximation has been used for the ma- 
trix formulation. Basis functions Ni,Ni , . . . ,N n have been chosen from a finite 
dimensional (of dimension n) subspace of H l . Hence, the test and trial functions 
can be finitely approximated as 


4> = CiNi(x, r) 


(13) 


4>\ = djNj(x, r) (14) 

where c,’s and d/s are suitable scalars 2 . Substitution of equations (64) and (65) 
in (63) results in the matrix formulation of the problem 


[LL 




rJx dx dx dr dr 


K(. 

*7 


— [ N{V(I> 0 • n rdT + A\ f iV,~ r • n rdT 
Jr, J r« it 2 


(15) 


F' 


where is the t, j entry of the stiffness matrix [K'\ and Fj is the i th entry in 
the load vector {jP'}. 


2 Since <f> is being suitably interpolated between the nodes, d/s here imply nodal values of the 
mean flow potential. 
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2. Inlet flow through fore exit plane only In a procedure similar to that of 
problem I, the weak form of problem II is 


Us 


50 2 50 50 2 50 


I* dx dx dr dr ” Jr f 

where P = r„ U T/ U r c and, 0 2 and 0 are from H 1 . 


) rdxdr = U fl f 0 rdT — A 2 f 0-j^-r • n rc (16) 

J V t JToo k 


As in problem I, a standard Galerkin finite element approximation has been 
used for the matrix formulation. The matrix formulation of problem II yields 

dj = U h [ Ni rdT -At f N~r . n rdT 

JTf JToo K 


[LI. 


t dNi dNj , dNidN ix , , 

(“sr ai" + sr-sr) rdxdr 


dx dx dr dr 


K ; 




( 17 ) 


where Kij is the i, j entry of the stiffness matrix [K] and Fi is the i th entry in the 
load vector {F}. 

Details of the stiffness matrix and load vector calculations are dealt with in 
the next sub-section. The constants A\ and A 2 for the problems I and II are 
found by imposing the compatibility condition which balances the flux across the 
boundaries. For problem I, it balances the flux across the nacelle, centerbody and 
fan face with the flux across Too, i.e., 


or, 


[ V0 O • nrdT = Ai f -r*nrdT 
J r» J r„ R 2 

/r* V0 O • nrdr 


A\ — 


( 18 ) 


/ r<so ^ 3 -r • n rdT 

For problem II, it balances the flux across the fan face with the flux across 
the outer boundary Too, i.e., 


U fl [ rdT = A 2 f 4rr»nrdT 
h { J r„ it 2 


or, 


Uh It, 

/rco^ r#nr<ir 


A,= 


( 19 ) 
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The constant -A s is calculated likewise. 

FINITE ELEMENT CALCULATIONS 

The global stiffness matrices and the global load vectors as defined in equa- 
tions (66) and (68), can be written as the composition of the element stiffness 
matrices and element load vectors respectively. Therefore, for example, in prob- 
lem II, we can write 

I*] = fm 

{f} = f ;{ f‘} 

where n e is the number of elements in the domain. The element stiffness matrix 
[K e \ and the element load vector {F*} for the problem are given by 


* - 


F; = U A / r< N! rdT - A, / r _ N/±r . n riV 


where / ne is the surface integral over the domain of the element, / r « and f r ^ are 
line integrals along element boundaries on the fan face and the outer boundary 
respectively, and N* is the shape function of the i th node of the element. 

Finite element calculations are done based on a parent element with local 
coordinates f and t] as shown in Figure 4. The element shape functions N f corre- 
sponding to each node i in the parent element are standard functions and therefore 
known. 


1. Surface Integrals To perform the finite element calculations on a parent 
element, an element map is constructed. The transformation under which each 
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element ft„ in the mesh is the image of a fixed parent element under a coordinate 
map T t is constructed as 


nodes 


T t :x— x iNl{S>r)) 

t=l 

(20a) 

nodes 

T e : r — £ r.-JV?(f,u) 

(206) 


*=1 


where nodes is the number of nodes on the parent element, and is 8 or 9 depending 
on whether it is an eight or nine-node element. The element H e to which T e maps 
the parent element is completely determined by specifying the x, r coordinates 
(x,-, r f ) of all the nodal points of f 2 e . Element shape functions iV*(x, r) are simply 
obtained from standard parent element shape functions N?(c,r)) by 

N-{x,r) = N-(i(x,r),r]{x,r)) (21) 


The derivatives of shape functions are obtained by the chain rule of differentiation, 


According to the element map, 
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By using the above relations (equations (20) through (23)), the element stiff- 
ness expression if,?- may be expressed as 

[ f{x,r) rdrdx = [ g(g, rj) r(g,rj) J(f, r))dgdr) (24) 

J n, J n p 
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where is the domain integral over the parent element and J(f, 77 ) is the Jacobian 
of the transformation T e given by 


. _ dx dr dx dr 

Ihidc 


( 25 ) 


A standard 4x4 Gaussian quadrature rule has been used to evaluate the integral. 
It is important to note that the mean flow calculations have been done both with 
and without the wave envelope elements. In one case, no distinction has been 
made between the elements in regions I and II and the elements in region III. All 
the regions have isoparametric rectangular finite elements. In another case, shape 
functions for the wave envelope elements in region III, differ from the rest in the 
mesh due to the fact that they simulate the inverse square decay behaviour as 
expected in a field due to a simple source. Mathematically a shape function may 
be expressed as 


Nr = n; ( 2 ) 


( 26 ) 


where N‘ is the standard shape function of node * at radius r,-. The results in 
both cases were virtually identical for the flow velocity we are concerned with. 
For compatibility with the radiation calculations the wave envelope elements were 
retained. 


2 . Boundary Integrals Three-node quadratic line elements lie along the 
boundaries and the generation of their topology and their nodal connectivity have 
already been discussed in the mesh generation scheme. The calculations of the 
line integrals are carried out by integrating along those sides of the parent element 
that are mapped onto the sides T e of the actual element fl e along which natural 
boundary conditions are prescribed. For definiteness, it has been assumed that the 
side f = 1 has been mapped onto T 8 . The line integrals have been parametrized 
with respect to rj. 
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The shape functions used for the line integrals are identical to the standard 
shape functions for the three-node Lagrangian line element. Element maps are 
created as discussed before and the elemental arc length is found by 

dr = Vdx 2 4- dr 2 


or, 

dv < 27 > 

' v ' 

j 

where j is the jacobian of the transfomation of 77 onto the arc length parameter 
in the x-r plane. The dot products, V<f> 0 • n for problem I and r • n for problems 
I through III need to be computed at each node on the relevant boundaries to 
evaluate the line integrals. Note that the constants A\ through A 3 need to be 
evaluated before constructing the load vector. 

D. THE SOLUTION PROCEDURE 

All of the boundary conditions are of the Neumann type and the differential 
equation is the Laplace equation. Hence, there is no unique solution to the prob- 
lems unless a reference value of the mean flow velocity is specified at any point 
in the domain. This does not affect the results because we are interested in the 
derivatives of the potential and not in the absolute values of the mean flow poten- 
tial. By penalization, the potential has been made zero at an arbritary point on 
the outer boundary Too. This penalization has been made at the elemental level. 
When the stiffness matrix of the element which occupies the penalized node is 
calculated, a very large value (1.0el5) is added to the diagonal entry in the matrix 
corresponding to that boundary node. Hence the velocity potential at that node 
is forced to zero after solution. The penalized stiffess matrix for that boundary 
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element looks like the following : 


K u • 

. • K lm ... 

K ln ‘ 

ifml 

• • K mrn + ~ « * • 

^ mn 

. K nl . 

. . K nrn 

1 

C 


where m is the penalized node number (local) and e (l.Oe-15) is the penalty param- 
eter. As a check, penalization was carried out at a different point in the domain, 
and the solution was found to differ from the previous solution by an arbitrary 
constant only. 

Since the penalization is carried out at the element level, the stiffess matrix 
and the load vector are never stored in assembled form. As each element stiffness 
matrix and load vector is formed, it is written down onto disk along with its nodal 
connectivity. The frontal solution method of Bruce Irons has been used to solve the 
algebraic system of equations [K]{(f>} — {F}. The principles of this technique are 
implied by the Gaussian process of forward elmination and back substitution. The 
frontal process alternates between accumulation of element coefficients (assembly) 
and elimination. Whenever an element is assembled its nodes are kept in active 
storage until their elimination. The active in-core storage at a point of time 
depends only on the “frontwidth” (number of active nodes at that time) which is 
much smaller than the dimension of the assembled matrix. This drastic reduction 
in in-core storage is the most important aspect of this scheme. Details of the 
scheme are, however, not discussed here. 

E. SUPERPOSITION OF THE SOLUTION FROM THE THREE PROBLEMS 

After the solutions to problems I through HI are obtained, they are added 
to the exterior field velocity potential to obtain the overall mean flow velocity 
potential of the flow field. Solution to the problem IV is the uniform flow field 
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whose velocity potential is given by 

<£u = U 0 X (28) 

The overall mean flow velocity potential is found out by pointwise addition of 
the solution vectors from the four flow problems 

<i>o — 4>i + 02 + <i>z + (29) 

However, the solutions to the problems I through IV have been obtained by an 
input of unit velocity at the exit planes, fore and aft, and in the exterior flow field 
i.e. 17/, = 1, U/ 3 = 1, and U Q = 1 respectively. Therefore, if the exit planes ( fore 
and aft) flow Mach number and exterior flow field Mach number are A//,, M/ a , 
and M u respectively, the superposed solution is found by 


4>o — + 4> u ) + + Mf 2 (j) 3 


(30) 



